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ABSTRACT 


The focus of this dissertation is on advanced development of the 
boundary element method for elastic and inelastic thermal stress 
analysis. New formulations for the treatment of body forces and 
nonlinear effects are derived. These formulations, which are based on 
particular integral theory, eliminate the need for volume integrals or 
extra surface integrals to account for these effects. The 
formulations are presented for axisymmetric, two- and three- 
dimensional analysis. Also in this dissertation, two-dimensional and 
axisymmetric formulations for elastic and inelastic, inhomogeneous 
stress analysis are introduced. The derivations account for 
inhomogeneities due to spatially dependent material parameters, and 
thermally induced inhomogeneities. 

The nonlinear formulations of the present work are based on a 
incremental initial stress approach. Two inelastic solutions 
algorithms are implemented: an iterative; and a variable stiffness 

type approach. The Von Mises yield criterion with variable hardening 
and the associated flow rule sire adopted in these algorithms. 

All formulations are implemented in a general purpose, multi- 
region computer code with the capability of local definition of 
boundary conditions. Quadratic, isoparametric shape functions are 
used to model the geometry and field variables of the boundary (and 
domain) of the problem. The multi-region implementation permits a 
body to be modeled in substructured parts; thus dramatically reducing 
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the cost of the analysis. Furthermore, it allows a body consisting of 
regions of different (homogeneous) material to be studied. 

To test the program, results obtained for simple test cases are 
checked against their analytical solutions. Thereafter, a range of 
problems of practical interest are analyzed. In addition to 
displacement and traction loads, problems with body forces due to 
self-weight, centrifugal, and thermal loads are considered. 
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A short list of notation is given below. All other symbols are 
defined when first introduced. A few symbols have different meaning 
in different context, but no confusion should arise. Matrices are 
indicated by bold print throughout this dissertation. 


A b ,B b ,C^ Boundary system matrices in assembled form 


k a ,B a ,C a Coefficient matrices of the stress equations in 

assembled form 
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E Young's modulus (in the Appendix, E represents the 

complete elliptic integral of the second kind) 
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Gravitational acceleration 


G ij’ F ij’ B ijk Kernels of the displacement equation 
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h The slope of the uniaxial equivalent stress, equivalent 
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Shear modulus 
Poisson's ratio 

Refers to coordinates of a field point 
Archimedes number 
Mass density 
(Real) stress 


xviii 



(Corrective) plastic part of initial stress 
Elastic stress 

Inhomogeneous part of initial stress 
Initial stress 

Thermal part of initial stress 

Yield stress * * 

System vector of initial stress 

t 

Fictitious particular integral density 

u Angular velocity 

Superscript 

c Components related to complementary function 

p Components related to particular integral (not to be 

confused with plastic component) 

. (Incremental) time rate of change 

Subscript 

, Spatial derivative 

i,j,k Indicial notation 

i,j,k = 1,2 in two-dimensions 
i,j,k = 1,2,3 in three-dimensions 
i,j,k = r,z in axisymmetry 

r,6,z Directions of cylindrical components 

xix 





CHAPTER 1 


INTRODUCTION 


1.1 GENERAL REMARKS 

1.2 THE HISTORICAL DEVELOPMENT OF BEM IN STRESS ANALYSIS 

1.3 SCOPE OF THE PRESENT WORK 


1 



CHAPTER ORE 


INTRODUCTION 


1.1 GENERAL REMARKS 

Through the years, scientists have formulated mathematical 
equations to describe the behavior of many physical phenomena in the 
field of continuum mechanics. However, due to the complexity of these 
equations, closed-form solutions are unobtainable for all but the 
simplest problems. In developing applications, engineers must 
therefore resort to approximate techniques for the solutions to 
physical problems. 

In the 1950s, the advent of the digital computer spurred the 
development of new approximate techniques or 'numerical methods'. 
Essentially, three numerical techniques evolved: the finite 

difference method; the finite element method (FEM); and the boundary 
element method (BEM). 

The finite difference method, which employed 'difference 
equations' for the generation of equation system, was the forerunner 
of these methods. However, the large size of the resulting equation 
system and its inability to readily deal with irregular boundaries 
limited its success. 

A more powerful and popular approach known as the finite element 
method (FEM) emerged. In this procedure the domain is divided into 
discrete elements, and trial functions are used to approximate the 
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functional variables across each element. Variational principles are 
applied to obtain a best fit solution for an appropriate set of 
boundary conditions. This method is very proficient in elastic and 
inelastic stress analyses and reasonable success has been achieved in 
dynamics. 

Owing to their mathematical and numerical simplicity, the 
development of these two methods was rapid. In contrast, the 
development was slower for a more complex technique known as the 
boundary element method (BEM). In recent years, however, researchers 
have focused considerable attention on this method mainly due to its 
advantages over the former methods. These advantages include: the 
ability to solve three-dimensional problems with greater efficiency, 
higher resolution results for stress concentration problems, increased 
accuracy, and greater ease in application for problems of infinite or 
semi-infinite regions. Moreover, the system equations are written 
only at nodal points on the surface of the body of interest, rather 
than throughout the entire domain. This leads to a set of equations 
that is smaller than that of competing methods. 

Two distinct, yet equivalent (Lamb, 1932) boundary element 
formulations exist. An indirect method, which introduces a ficticious 
set of functions as an intermediate step in the solution process, and 
a more popular direct formulation, which utilizes all real, physical 
variables. Essentially, both methods consist of transforming the 
governing differential equation into a boundary integral equation. 
The surface of a body is divided into boundary elements, and shape 
functions are used to represent variation of quantities across the 
elements in terms of their nodal values. These integral equations are 
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integrated numerically, generating a system of equations at boundary 
nodes. Standard solution procedures are employed to obtain results 
for a prescribed set of boundary conditions. 

1.2 THE HISTORICAL DEVELOPMENT OF BEM IN STRESS ANALYSIS 

The first rigorous work on integral equations was published by 
Fredholm in 1903. Since that time, extensive research has been 
carried out by a number of researchers such as Kellog (1929), 
Muskhelishvili (1953), Kupradze (1964) and Smirhovi (1964). The most 
notable work related to elastostatics was conducted by Mikhlin (1957, 
1964,1965). However, due to the complexity of finding analytical 
solutions, most of these early works were of mathematical nature 
dealing in topics of existence and uniqueness. 

The advent of the digital computer brought about a change in the 
usefulness of these mathematical formulations. The first application 
in elastostatics by Jaswon and Ponter (1963) dealt with torsion in 
elastic bars. A general elastostatic (direct BEM) approach based on 
the displacement formulation derived from the Somigliana's identity 
(1885) was established by Rizzo (1967). Further contributions were 
made by Cruse (1969, 1973), and Lachat and Watson (1976). An 
alternative, indirect, elastostatic formulation was developed by 
Banerjee (1969), Butterfield and Banerjee (1971), Watson (1973), 
Tomlin (1973), and others. 

As the method developed, the analysis was extended to other 
areas. Cruse (1967) pioneered the study on transient elastodynamics, 
and Rizzo and Shippy (1977) laid forth an efficient procedure to 
account for steady-state body forces such as centrifugal and thermal 
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loadings. Kermanidis (1975) produced the first elastic, axisymmetrie 
formulation, followed by Cruse, Snow and Wilson (1977) who extended 
the axisymmetrie analysis to include thermal and centrifugal body 
force effects. And in 1979, Nigam (1979) introduced an axisymmetrie 
formulation capable of handling non-axisymmetric boundary conditions. 

It soon became apparent that material nonlinearities, in the form 
of initial stress (or initial strain), could be incorporated into the 
elastic analysis through a volume integral in a manner analogous to 
body force, and therefore, quasi-static algorithms for plasticity and 
creep analysis were developed. Swedlow and Cruse (1971) presented the 
first elastoplastic BEM formulation, followed by Riccardella (1973) 
who is credited with the first two-dimensional formulation. Further 
development by other workers followed: Rzasnik and Mendelson (1975), 
Mendelson and Albers (1975), and Chaudonneret (1977). These authors, 
however, overlooked the strong singularity present in the domain 
integral of the interior stress rate equation. In 1978, Bui (1978) 
presented a correct treatment of this integral and indicated the 
existence of the free term. Nevertheless, this integral was strongly 
singular and the numerical integration was still difficult. 

Banerjee and co-workers developed a two-dimensional (Banerjee and 
Mustoe, 1978), elastoplastic BEM formulation and were the first to 
extend the analysis to three-dimensional (Banerjee, Cathie, and 
Davies, 1979) and axisymmetrie (Cathie and Banerjee, 1980) media. And 
in 1982, Cathie and Banerjee (1982) presented a time independent 
inelastic analysis. Their work is unique in the sense that the stress 
rates are calculated via numerical differentiated displacement rates, 
and therefore, the strongly singular domain integral is avoided. The 
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method is computationally efficient, however, some accuracy is lost in 
the numerical differentiation. 

Mukherjee and co-workers also have concentrated considerable 
effort in this area. Their contributions included: a time-dependent 
inelastic analysis using power creep; and the constitutive relations 
due to Hart (Mukherjee and Kunar, 1978); and the implementation of an 
axisymmetric visco-plastic analysis (Sarihan and Mukherjee, 1982). In 
addition, they have presented an integral equation for stress rates 

4 * 

written as a function of a gradient of initial strain rates which 
eliminates the appearance of the strongly singular Lebesque integral 
over the domain in favor of a weaker volume integral. 

Other researchers, such as Telles and Brebbia (1981), Kobayashi 
and Nishimura (1980), and Telles (1983), have presented formulations- 
and have solved a, variety of plasticity problems. 

The most notable advancement in stress analysis via the boundary 
element method is the development of the general purpose, three- 
dimensional, dynamic and inelastic, stress analysis program - BEST3D - 
developed for the National Aeronautics and Space Administration (NASA) 
by Banerjee, Wilson and Miller (1985), Banerjee and Ahmad (1985) and 
Banerjee and Raveendra (1986). Highlighted in this work is the use of 
isoparametric quadratic shape functions for modeling the yariation of 
functions over the boundary elements and volume cells. In addition, 
the system includes a sophisticated numerical integration scheme 
(Banerjee, Wilson, Miller, 1985) and a multi-region facility. 

The BEST3D system contains the most advanced plasticity analysis 
to date. A Von Mises model, a two-surface model and a thermally 
sensitive anisotropic plasticity model, are all contained in the 
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system. Furthermore, an iteration acceleration scheme is used to 
reduce the number of iterations needed for convergence. 

All the aforementioned plasticity formulations are based on 
iterative procedures that work successfully, but often takes unduly 
large number of iterations to converge to the correct solution, 
particularly in problems involving a high degree of nonlinearity such 
as the loading close to the collapse state of stress when significant 
amount of plastic zones develop. It was to this end that Banerjee and 
Raveendra (1987) presented the first direct or 'non-iterative' two- 
dimensional elastoplastic analysis which is comparable to the variable 
stiffness method in the finite element analysis. 

In any thermal stress analysis, it is important to account for 
the variations in elastic modulus. In 1968, Rizzo and Shippy (1968) 
presented a BEM formulation for inhomogeneous elastic inclusions. 
More recently Tanaka and Tanaka (1980) introduced a thermoelastic 
formulation for inhomogeneous material. In this analysis the material 
parameters are independent of temperature, and vary only as a function 
of position. Finally, Ghosh and Mukherjee (1984) presented a two- 
dimensional iterative formulation for thermoelastic deformation of 
inhomogeneous media. In their work, they assumed the shear modulus 
varied linearly with temperature. In spite of these efforts the 
problem of elastic inhomogeneity remains essentially unresolved. 

1.3 SCOPE OF THE PRESENT WORK 

Once considered an analysis for mathematicians, the boundary 
element method is now accepted in the engineering community as a 
powerful and versatile tool for solving practical problems. In many 
applications, the boundary element method has proven itself superior 
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to other numerical methods, particularly in three-dimensional linear 
analysis. However, further refinement is necessary if inelastic 
boundary element analysis is to assert a similar claim. The primary 
objective of the present work is to rectify some of the deficiencies 
in this area. Furthermore, in this dissertation, considerable effort 
is put forth in extending all formulations to their axisymmetric form. 
The realization that many stress analysis problems of industrial 
interest are axisymmetric in nature, and the fact that the 
axisymmetric analysis costs a fraction of the three-dimensional 
counterpart, makes the axisymmetric analysis an invaluable tool. 

To better appreciate the new boundary element formulations 
presented in this dissertation, the conventional BEM formulations for 
elastic and inelastic thermal analyses are first presented in Chapter 
2. Included is a comprehensive look at the details of axisymmetric 
inelastic stress analysis. 

The need for only surface discretization is a significant 
advantage BEM has over other methods requiring full domain 
discretization. However, this advantage is partially diminished in 
thermal and inelastic analysis where volume integration is required. 
To rectify the matter, a new formulation, based on particular integral 
theory is developed. This new procedure, incorporates the body force 
effect into the boundary element system without additional surface or 
volume integration. Application of the method for gravitational, 
centrifugal and thermal body forces is presented in Chapter 3 and in 
Chapter 4, the formulation is also extended to the inelastic analysis. 

In Chapter 5 one of BEM's greatest deficiencies is addressed; 
that being its inability to deal with material inhomogeneities. Both 
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natural material inhomogeneities as well as thermally induced 
inhomogeneities are considered for both linear and non-linear thermal 
analysis. 

All formulations of the present work are implemented into a 
general purpose, multi-region system capable of local definition of 
boundary conditions. Quadratic isoparametric shape functions are used 
to model the geometry and the field variables of both the boundary 
element and volume cells. In Chapter 6, a variety of linear and 
nonlinear problems are presented. In addition to displacement and 
traction boundary loads, body forces due to self-weight, centrifugal, 
and thermal loads are considered. In most problems, the body of 
interest is sub-structured for a multi-region analysis. This 
dramatically reduces the time and cost of the analysis as will be 
discussed. 
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CHAPTER TWO 


CONVENTIONAL BOUNDARY ELEMENT FORMULATION 
FOR ELASTIC AND INELASTIC STRESS ANALYSIS 

2.1 INTRODUCTION 

In this chapter, the conventional formulation of the boundary 
integral equations for elastic and inelastic thermal stress analysis 
is presented. In addition to the two- and three-dimensional 
formulations, the axisymmetric case is also developed. The body 
forces and nonlinear effects are incorporated in these formulations in 
the conventional manner through volume integrals or addition surface 
integrals. 

After a brief discussion on the elastoplastic constitutive 
relations, two BEM algorithms for the solution of inelastic problems 
are described. The first algorithm is an iterative procedure which 
includes a time saving feature which reduces the number of iterations 
needed for convergence by utilizing the past history of initial stress 
rates to estimate the values of the initial stress rates of the next 
load increment. The second algorithm is a direct procedure similar to 
the variable stiffness approach in the finite element method. This 
procedure exploits certain features of the constitutive relationship 
to express the unknown nonlinear initial stress rate tensor as a 
scalar quantity which then can be eliminated from the boundary 
equation system through a back substitution of the (modified) stress 
rate equations. 
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The axisymmetric, two- and three-dimensional formulations are 
implemented in a multi-region code which utilizes quadratic 
isoparametric shape functions to model the geometry and field 
variables of the boundary and domain of the body. Numerical 
implementation procedures are briefly discussed, including techniques 
that are used to calculate coefficients of the integrals over strong 
singularity points on the boundary and domain of the body. Finally, a 
number of examples are included to demonstrate the accuracy and 
convergence of the plasticity algorithms. 


2.2 BOUNDARY INTEGRAL FORMULATION 
2.2.1 Governing Equations 

The governing equation of stress for a body in equilibrium with- 


domain V and surface S can be expressed in a cartesian system x A as 


a . . . + f . a o' (2.1) 

ij» J i 

i,j = 1,2 for two-dimensional 
i,j = 1,2,3 for three-dimensional 
where a.. = .(x) is the stress tensor, 

1 J 1 J 

f. = f.(x) represents the total contribution of all body 

forces present, and 

the comma indicates spatial derivatives. 


Throughout this dissertation, when the dimension of the problem is 
irrelevant to the point of the discussion, three-dimension will be 
assumed. Solutions of equations (2.1) are subject to boundary 
conditions: 
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u i = U i (x) 

on Sj 

(2.2a) 

t i = ff ij n j = T i (x) 

on S 2 

(2.2b) 

a compatible combination of the two 

on S 3 

(2.2c) 


where S = + Sj + S^, 

is the displacement vector, 
t i is the traction vector, 

n j is the boundary normal vector, 

U^x) are the imposed displacements, and 

T^(x) are the imposed tractions. 


In equation (2.2b) we have used the Cauchy traction relation 


t. 

1 


T ij n j 


(2.3) 


and displacements are related to strain via 


e . . 
1J 


1/2 (u i.j +u j.i ) 


(2.4) 


The (total) strain tensor e. . can be decomposed into an elastic 

J 

strain e®^ an d an initial strain 

8 ij = 8 ij + 8 ij (2.5) 

The elastic strain is directly related to stress through the elastic 
constitutive relations 


0 iJ " D ijkl 8 kl 


(2.6a) 


e 


e 

ij 


= C 


e 

ijkl ff kl 


(2.6b) 


where 


D ijkl “ 6 ik s jl + x 6 ij 8 kl 


(2.7a) 
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(2.7b) 


0 1 + V V 

C ijkl = IT 5 ik 8 jl ' E 8 ij 8 kl 


is the Kronecker delta 
X and n are Lame constants 
E is the modulus of elasticity, and 
v is Poisson's ratio 


Note: D ®j kl c kl pq " 8 ip 8 jq 


( 2 . 8 ) 


D? and C® .. are symmetric with respect to (ij) and (kl) 

lJKl lJKX 

The initial strain, by definition, is the difference between the 
total strain and the elastic strain. The initial strain occurs from 
various effects, such as plastic deformation, thermal loading, or 
strain present in a body before loading occurs. In a thermoplastic 
analysis, initial strain will be defined as the summation of plastic 
strain e?. and thermal strain 

X J X J 


o° P , t 

e ij - e ij + 8 ij 


(2.9) 


Substituting equation (2.5) into equation (2.6) yields a relation 
between total strain and stress 


a ij ~ D ijkl 8 kl " 


(2.10a) 


P e . o 

c ijkl ^kl + 8 ij 


(2.10b) 


where 


TN e 0 

D ijkl e kl 


(2.10c) 


Substituting equation (2.4) and (2.10a) into equation (2.1) 
produces the governing differential equation of equilibrium in terms 
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of displacement, body force, and initial stress: 


( X*m> Uj>lj + „ u UJ) * f ± = »2 Jj3 (2.11) 

This is the Navier equation with an initial stress term present. 

2.2.2 Two- and Three-Dimensional Integral Formulation 

In deriving the boundary integral equations for general elastic 

and inelastic analysis, we consider an elastic body under two distinct 

equilibrium states. The first state (u^ f A , o^) will be 

considered the real state, and the second state (u., t., f-.e?-. a- ■) 

1’ x’ x* Xj' °Xj' 

will be considered an arbitrary state. 

The following integral statement can be inferred from equations 
(2.6) and symmetry property of 

J •« 'ij dv - 1 «u dV (2 - 12 > 

V V 

Using equation (2.5) this equation can be rewritten in terms of total 
strain as 


J °ij 8 ij dV - J e ij ®ij 

V V 


D 


e 

ijkl 


Using equations (2.6b), 
we can show 


dV + J ejj dV (2.13) 

V 

(2.10c), and the symmetry property of 


J «1J E ij dv - J H j 4} dT < 2 - 14 » 

V V 

Substituting this result in equation (2.13) gives 
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(2.15) 


1 «u av = J 4j «ij dV * S ‘h °h dV 

v v V 

Employing equation (2.4) and utilizing the symmetry property of 

the stress tensor, equation (2.15) is transformed into 


I Hi "i,J dv - I "i.j «ij + I Hi Hi dv 

V v V 


(2.15) 


Applying the divergence theorem to the first two terms of the 
above equation and substituting equations (2.1) and (2.3) into this 
result we arrive at 


J t i (x) u ± (x) dS + j f ± (x) u^x) dV = j u|(x) t i (x) dS 

S V S 

+ J uj(x) f^x) dV + j ejj dV (2.17) 

V V 

In the absence of initial stress, this equation reduces to the 

celebrated Betti's reciprocal work theorem. 

We now define the arbitrary, elastic stress state to be 
equivalent to the state of stress given by the Kelvin point force 
solution for point force e k (£) acting at ? i . 


^(x) = 6(x,£) 6 ik e k (£) 

u i (x) " G ik (x,£) e k ($) 

«®j(x) = B iJk (x,5) e k ($) 

t i( x) = F ik (x£) e k ( ^ ) 

where 

8(x,£) is the Dirac delta function, and 
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G ik» B ijk» and F ik represent the displacement, strain and stress 
at point x^^ i n an infinite space due to a unit point force at 
(both the two- and three-dimensional functions are defined in the 
appendix). Noting the properties of the Dirac delta function, we 
can show that 

J f^x) u i (x) dV = J 8(x,lj) 8 ik e k (5) u^x) dV = u k (£) e k (J;) 

(2.19) 

Substituting equations (2.18) and (2.19) into equation (2.17) and 
using the orthogonality property to separate the three independent 
components (and cancelling the e k (£)), the displacement at point 
can be expressed as 

C ij(£) = J [G i j(x,^) t i (x) - F i j(x,?) u i (x)] dS(x) 

S 

+ | G i j(x,5) f ± (x) dV(x) + J* B ik j(x,£) o° k (x) dV(x) 

V V (2.20) 

where C^j(£) = 8 .jj for a point in the interior of the domain. When 
point ^ i S brought to the boundary x Q , C^j(x Q ) must be derived from 
the singular treatment of surface integral involving the kernel. 
The resulting tensor function C is depended on the subtended angle 
of the tangent plane at x Q . For a point on the smooth surface = 

8^. The surface integral involving kernel F^j in equation (2.20) 

must be treated as a Lebesgue integral. 

The integral equation for the strain at an interior point is 

found analytically by substituting equation (2.20) (with C. . = 5..) 

1J 1J 

into the strain-displacement relations, where differentiation is with 
respect to the field point % : 
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e.jU) = J [G^ ij (x,5)t k (x)-F^ ij (x,5)u k (x)]dS(x) 

S 

+ j G^ij(x,5)f k (x)dV(x) + J B kli j(x,§)<y kl (x)dV(x) (2.21) 

V V 

By introducing this result into the stress-strain equation (2.10a) the 
stress integral equation is derived: 

0^(5) = J [Gj ij (x. 5 )t k (x)-Fj ij (x^)u k (x)]dS(x) 

S 

* * 

+ J* G^ i j(x,^)f k (x)dA(x) + | B^ li j(x,^)«T kl (x)dV(x) (2.22) 

V V 

The kernels , F^j and are defined in the Appendix. 

The last integral in equations (2.21) and (2.22) is strongly singular 

and must be treated as Lebesgue integral. Treatment of these 
* 

integrals is discussed in subsequent sections. 

Due to strong singularities in the kernel functions, equations 
(2.21) and (2.22) are not used for calculations of stresses and 
strains on the boundary. Instead, an alternate procedure will be 
presented in later sections. 

2.2.3 Axisymmetric Integral Formulation 

The axisymmetric boundary integral equations can ,be derived 
from the three-dimensional BEM equation (2.20) which satisfies the 
governing differential equation (2.11). In principle this is 
accomplished by: recasting the three-dimensional BEM equation in 
cylindrical coordinates; applying appropriate tensor transformations; 
and integrating analytically to remove the angular (6) dependency. 
However, once the axisymmetric displacement kernel G^ is known, the 
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F^j and kernels ean be derived an alternate way starting with the 

solution (the G_ kernel) for the displacements due to a ring load 
intensity (Cathie and Banerjee, 1980). The strain due a ring load 


(the B_ k kernel) is obtained by substitution of the kernel into 
the axisymmetric strain-displacement equations, where differentiation 
is with respect to x. The surface traction due to a ring load (the 
F ij kernel) is found by introducing the B ijk kernel into Hooke's law 
and multiplying by appropriate normals. 

The axisymmetric form of the displacement integral equation 


can then be expressed as: 


^ij ( V 


= J [G i j(x,^)t i (x)-F i j(x,^)u i (x)]dC(x) 
C 


+ J G i j(x,^)f i (x)dA(x) 
A 



B ik j (x,£) cr? k (x)dA(x) 
i,j,k = r,z , 


(2.23) 


where 


u ^, t^, fj_, and <Tj_ k are expressed in radial and axial components. 


^ij ~ ^ij f° r i nter i° r points and is dependent on the surface 
geometry at \ for boundary points, and 


G 


ij 



and B ijk 


are defined in the Appendix. 


Note that the surface integration is over a curve dC and the 
volume integration is over an area dA. Furthermore, the 2nR term 
which appears after integration in the 6 direction has been absorbed 
in the kernels. The domain integral need only be evaluated where the 
initial stresses are non-zero. 
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Axi symmetric integral equations for stress and strain at interior 
points are derived in a manner analogous to the two— and three- 
dimensional case. 

2 .3 NUMERICAL IMPLEMENTATION 

The integral equations discussed in the preceding sections become 
of practical interest only when numerical techniques are employed for 
their solution. It should be noted that these integral equations are 
exact statements, and any error in the final BEM solution is caused by 
approximations in the numerical implementation. 

The numerical implementations of the integral equations consist 
of: discretization of the surface of the body (and domain when 
necessary) and discretization of the integral equations; numerical 
integration; assembly, and solution of these equations. 

2.3.1 Discretization and Numerical Integration 

Discretization - The body under consideration is divided into 
discrete boundary elements and, when applicable, volume cells. 
Quadratic isoparametric curvilinear shape functions, presented in 
Appendix I, are used to approximate the geometry and the field 
variables over the boundary elements and volume cells in terms of 
their nodal values. After discretization, the boundary displacement 
equation can be expressed in the following manner. 

CiiVs) = 5 [ I Gij(x,5)N Y (n)dS m ] t^ 
m=l s m 
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(2.24) 


P 

+ H 1 Gii(x,5)N P (ii) dvP] fjPP 

P-1 V p 

P 

+ }[ I B ikj (x.?)N P (n)dvP] (<r° k )pp 
P-1 V P 

where i,j,k = r,z 

M = number of boundary elements, 

P = number of volume cells, 

N y (ti) represents shape function of boundary elements, and 

nP(ti) represents shape function of volume cells. 

Summation over y and 0 is implied. 

The bars indicate nodal values and the integration coordinate x 
has been expressed in local coordinates -n via the above shape 
functions. 

The stress equation can be discretized in a similar fashion: 

V {) - H J oJ 1J (x.?)B T (n) 

m-l s n 

“HI Fj ij <x,?)N T (n) dS“] u l m 
m=l s “ 

+ 2 [ J oSij(x.5)NP(,) dVP] 

P=1 V P 
p 

+ 5 [ J Bj lij (x,5)N P (r l ) dV p ] (<y kl ) Pp (2.25) 

P-1 v p 

A similar discretization of the strain equation is possible, 
however, the implementation is unnecessary. Once the stress state at 
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a point is known, the strains can be calculated by a direct 
application of the constitutive relation. 

Numerical Integration - The complexity of the integral in the 
discretize equation necessitates the use of numerical integration for 
their evaluation. The steps in the integration process for a given 
element is outlined below: 

1. Using appropriate Jacobian transformations (given in 
Appendix I), the curvilinear line element, surface element 
or volume element is mapped on to a unit line, a plane unit 
cell or a three-dimensional unit cell, respectively. 

2. Depending on the proximity between the field point (£^) and 
the element under consideration, there may be element sub- 
division and additional mapping for improved accuracy 
(Banerjee and Raveendra, 1986). 

3. Gaussian quadrature (Appendix I) formulas are employed for 
the evaluation of the discretized integral over each element 
(or sub-element). These formulas approximate the integral 
as a sum of weighted function values at designated points. 
The error in the approximation is dependent on the order of 
the (Gauss) points employed in the formula. To minimize 
error while at the same time maintaining computational 
efficiency, optimization schemes are used to choose the best 
number of points for a particular field point and element 
(Watson, 1979). 
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2.3.2 Treatment of the Singular Integrals 

The integration of the kernel functions are singular and require 
special attention when an integration point (x) coincides with the 
field point (£). 

The Displacement Equation - The integrations of the and 

kernels of the displacement equation are weakly singular under this 

circumstance and can be integrated numerically using element 

subdivisions near the singular point (Banerjee and Raveendra, 1986). 

On the other hand, the integration of the kernel of the 

displacement equation is strongly singular and must be integrated as a 

Lebesgue integral over the singular element. This integral can be 

decomposed into a free term and a Cauchy principal-value integral, but 

accurate numerical integration is still difficult (Banerjee and 

Raveendra, 1986). When constructing the equations for the boundary 

system, the singular integration involving the F kernel can be 

**■ 

circumvented by using the 'rigid body' displacement technique (Swedlow 
and Cruse, 1971) or the analogous 'inf lation mode, ' for axisymmetry 
(Nigam, 1979). A brief description of these procedures is given 
below. 

Rigid Body: For every boundary node of the system there is a 2 

by 2 (or 3 by 3 for three-dimension) block of coefficients on the 
diagonal of the F matrix, corresponding to the singular node, which is 
difficult to determine by numerical integration. Each of these terms, 
however, can be determined independently by assuming two (or three) 
different admissible displacement fields corresponding to a rigid body 
displacement in each of the two (or three) directions. All tractions 
(and body forces) are zero for this calculation. 
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Inflation Mode: Since a rigid body displacement is not 

admissible for the radial direction in axisymmetry, the two unknown 

coefficients (corresponding to the radial displacement) are determined 

by an inflation mode (i.e. a linear displacement in the radial 

direction), u = r and u„ = 0 . With body forces set to zero, the 
r z 

tractions related to this displacement are t = 2(X + fi) n r and t 2 
= 2Xn , in which n„ and n_ are normals on the boundary and X and 
H are Lam6 constants. The remaining two unknown coefficients in the 

l % 

displacement equation for the z-direction. are determined by using a 
rigid body displacement in this direction. 

It should be noted that when displacements are required on the 
boundary between nodal points, the values should be calculated via 
the shape function of the boundary element . 

The Stress Equation - The stress (or strain) calculation is 
handled in two ways depending on whether the field point falls on the 
boundary or in the Interior of the domain. 

Stress on the Boundary: For a point on the boundary, all kernel 

functions of the stress equation are all strongly singular and 
difficult to integrate numerically. However, the stress on the 
boundary can be obtained from the boundary tractions and 
discplacements without any integration as originally proposed by 
Cruse (1974). A slight modification must be made to incorporate the 
initial stress. In this procedure the stress-strain relations, the 
Cauchy traction equation and the equations relating local and global 
gradients of displacements are utilized in writing an expression for 
the boundary stress. Demonstrated for the axisymmetric case, these 
equations can be written in matrix form as follows 
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(2.26) 

where = \ + 2|i c 2 = x 

n and n are normals on the boundary, and 
r z 

u r#s and u z s are the local displacement gradients on the 
boundary. 

Using shape functions and their first derivatives the vectors on the 
right hand side of equations (2.26) can be expanded as a matrix 
multiplied by a vector consisting of nodal values of displacement, 
traction and initial stress. Equation (2.26)is then pre-mu 1 tip lied by 
the inverse of the matrix on the left hand side. 

Stress in the Interior: When the field point lies in the 

interior of the body, the and kernels of the stress 

equations are well behaved and require no special attention. However, 

the kernel is strongly singular and must be integrated as a 

Lebesgue integral over the singular point. Once again, this type of 

integral can be decomposed into a free term (given in the Appendix) 
and a Cauchy principal-value integral. Although difficult and 
expensive, reasonable accuracy can be obtained when this Integral is 
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numerically integrated using techniques described by Banerjee and 
Raveendra (1986). However, for high accuracy, the coefficients of the 
B® kernel corresponding to a singular nodal point should be 

3. J K X 

evaluated using a procedure known as the 'Initial Stress Expansion 
Technique, ' described below. 

2.3.3 Initial Stress Expansion Technique 

In this procedure, the coefficients of the stress equations 
related to the non-singular nodes are integrated in the usual manner. 
In each stress equation there remains three undetermined coefficients 
in 2-D (six in 3-D, or four in a axisymmetry) corresponding to the 
initial stress at the singular node. In a manner analogous to the 
'rigid body' technique, each of these coefficients are calculated by 
assuming one of the three (six or four) admissible stress states and 
compatible displacement fields given in Table 2.1 (2.2 or 2.3). For 
each stress state, one unknown coefficient in each stress equation can 
be determined. It should be noted that in order to apply this method 
the entire region must be covered with cells. In an iso-thermal 
analysis in which the plastic yield zone is small, this technique 
appears inefficient since it requires the presence of volume cells 
throughout the elastic region which otherwise would be unnecessary. 
Although this is true in a single region program, such is not the case 
in multi-region code since the technique is applied to each region 
independently and only used in regions where volume cells exist. 
Therefore, the zones in which plastic yielding is expected to occur 
are isolated in separate regions fully populated with cells, and the 
elastic regions remain free of any volume cells. 
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2.3.4 Singularities at the Origin in Axisymmetry 

Additional singularities exist in the axisynunetric kernel when 
the field point falls on the origin ($ = o). This problem is 

resolved by moving the field point a small radial distance from the 
origin. This is a common practice also employed in the finite element 
analysis of axisymmetric problems. In the equations for stress and 
strain, the minimum distance should be 0.5% of the difference between 
the maximum and minimum z coordinates at the origin. For the 
displacement equation this distance can be considerably less. 
Alternatively, a separate limiting form of the displacement kernels 
for a field point on the origin can be derived (Bakr and Fenner, 
1983) . 

Finally, for circumferential strain on the boundary at the 
origin, the relation = e r = 3u r /dr proves useful instead of the 
usual e e = Up/r. 

2.3.5 Assembly and Solution of Equations 

A displacement equation is written for each boundary node and a 
stress equation is written at points of interest. The resulting 
coefficients are modified appropriately for cases where the functions 
(boundary conditions) are referred to local boundaries. The 
displacement and stress equations are assembled by collecting the 
known and unknown values of tractions and displacements and their 
coefficients together. At the common interface of substructured 
regions, the equilibrium and compatibility conditions are invoked in a 
manner described by Banerjee and Butterfield (1981). The final system 
equations can be cast as: 
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(2.27a) 


A b x = B b y + C b o° 
a = A ff x + B^y + C®cr° (2.27b) 

where 

x is the vector of unknown variables at boundary and 

interface nodes. 

y is the vector of known variables, 

a° is the vector of initial stress (includes thermal load 

contribution) 

A b ,B b ,C b are the coefficient matrices of the boundary 
(displacement) system, and 

A ff ,B a ,C a are the coefficient matrices of the stress equations. 

It should be noted that A b is a square matrix. Furthermore, in a 
substructured system the matrices A b and B b are block banded 
whereas matrices C b , A*, B v and C a are block diagonal. 

Solution: Standard numerical procedures are used to solve for 

the knowns in equations (2.27a), after which the stresses in equation 
(2.27b) can be determined. 

The solver, employed in this work, is part of a software package 
from LINPACK (Dongarra, et al, 1979), which decomposes the A matrix 
into an upper triangular submatrix by the Gaussian reduction process. 
This procedure allows for an efficient resolution for vectors of 
different values on the right hand side of the same equation. This is 
essential for an efficient plasticity algorithm. 
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2.4 CONVENTIONAL BEM FORMULATION WITH BODY FORCES 

The general form of the governing differential equation for the 
deformation of a homogeneous isotropic body subjected to 
gravitational, centrifugal and thermal loading is given in equation 
(2.11) and is repeated here. 


(x+|i) 

U j.ji + * U i.jj 

+ f i " *ij,j 

(2.28) 

which 




f . 
i 

= fS + f? 

l 



f? 

l 

= P g e 3 

for gravity loading in z direction 

f? 

l 

= pw^(x2e2+X2e2^ 

for centrifugal loading about 

z axis 

a?. 

ij 

" f* T 8 ij 

for thermal loading 



p is the material density, g is the gravitational acceleration, 
is a unit vector in the i direction, w is the angular 
velocity, and x ^ is the spatial coordinate from the center of 
rotation, p = a (3>.+2fi), a is thermal coefficient of expansion, 
and T is the change in temperature. 

For simplicity, rotation is assumed to be centered at the origin about 
the z axis. The above equation is subject to boundary conditions 
given in equation (2.2). 

The three-dimensional boundary integral equation for displacement 
satisfying equation (2.28) is given in equation (2.20) and is repeated 
here. 

C ij($) u i(5) = J [G i j(x,§)t i (x) - F^x.SJu^x)] dS(x) 

S 

+ J G ij (x,5)f i (x) dV(x) + { B ikj (x,5)<xJ k (x) dV(x) (2.29) 

V V 
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and the corresponding equation for stress is given in equation (2.22). 

Two drawbacks exist in the formulation. First, numerical volume 
integration is time consuming and expensive, and second, the volume 
integrals of the stress equation are strongly singular leading to 
difficulties in the numerical integration. 

For these reasons, Rizzo and Shippy (1977) (and Cruse, Snow and 
Wilson, 1977 in axisymmetry) transformed the volume integrals of 

equation (2.29) into a surface integral, through an application of the 

* * 

divergence theorem (assuming a steady- state temperature distribution). 
The three-dimensional version is expressed as 

C ij ( ^ )u i (?) = I [G ij (x ’^ {t i (x) “ h(x) “ F i j(x,^)u i (x)] dS 

S 

dr d 

+ u f (h(x)+pT(x)} ~ — (h(x)+pT(x)} r .+3pa> n ,(x)r] dS 

o j on on 9 J 

S (2.30) 


where 

r is the eucidean distance between x and £, 

T(x) is the temperature change at x, 

n^(x) is the boundary normals at x, 

P o = (l-2p)/16nn(l-v) ; vis Poisson's ratio, and 
f ± = h, i (x) . 

The method is effective in many applications, particularly for 
gravitational and centrifugal analysis, but is restricted to steady- 
state temperature distribution and still requires additional surface 
integration. However, Rizzo (1977) and Shippy stated that the concept 
is applicable to diffusive systems and Masinda (1984) has, without 
evidence of implementations, produced formulations. Nevertheless, the 
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formulation is restricted to problems in which no heat sources are 
present and is only valid for uniform initial temperature 
distribution. For a general transient thermoelastic analysis, one 
still requires the use of volume internals. 

2.5 INELASTIC BEM FORMULATION BASED ON VOLUME INTEGRALS 

2.5.1 Incremental Theory of Plasticity 

Before proceeding with the inelastic BEM formulation, it is 
necessary to explore the constitutive relations that govern the 
material nonlinearities. The stress-strain relationships based on the 
'path independence in the small' theory of plasticity essentially have 
three components. 

1. A yield criterion - defines the limit of elastic behavior. 

2. A plastic flow rule - relates the irrecoverable plastic 
strain increment to the state of stress in a material. 

3. A hardening rule - defines the expansion or contraction of 
the subsequent yield surfaces resulting from continuous 
plastic flow. 

The nature of the yield surface and hardening parameters are 
material dependent. In this dissertation the Von Mises yield 
criterion with strain hardening, Drucker's postulate for associated 
plastic flow, and isotropic hardening is assumed. These assumptions 
provide a good model for ductile material, such as metal, under 
monotonic loading. 

The Von Mises Yield Criterion - The Von Mises Yield Criterion 
with strain hardening is assumed in the present work. This can be 
expressed in terms of the deviatoric stress tensor and plastic 
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strain e?. 

ij 


as 



<r Q = o 0 (e?j) is the uniaxial yield stress, and 

h = the slope of the uniaxial equivalent stress plastic strain 
curve . 


Plastic Flow Rule - The plastic flow rule (Drucker, 1952) states 
that the plastic strain rate tensor is linearly related to the 
gradient of the plastic potential Q through the stress point. This 
implies that the plastic strain rate tensor is normal to the surface 
of the plastic potential. This condition known as normality principle 
is generally expressed as 



(2.32) 


\ is non-negative scalar variable which depends on the current 
stress rates and the past history of loading. 


For metals, 'associated flow' (i.e., Q = F) is usually assumed 
and will be adopted in this work. Other materials, such as soil, 
require more complicated functions, and therefore 'non-associated 
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flow' (Q ^ F) is often assumed. 


Isotropic Hardening - A hardening rule defines the subsequent 
yield surfaces resulting from continuous plastic deformation. 
Different hardening rules have been proposed to model the behavior of 
various materials and loadings. 

The isotropic hardening rule which is employed in the present 
work, assumes the yield surface uniformly expands about the origin in 
the stress space while maintaining its shape, center and orientation 
during plastic flow. Figure 2.1 illustrates this behavior. During 
loading along path OA, the stress state is elastic. At point A, the 
elastic limit is reached and additional loading along path AB expands 
the yield surface to point B. Upon unloading, the behavior is elastic 
along path BC until point C is reached, beyond which yielding occurs 
only once. 

Let us assume a yield function for isotropic hardening material 
can be expressed as 


F ejj, h) = 0 (2.33) 

where cr^ is the current state of stress, is the total plastic 
strain and h is the hardening that may vary with plastic strain. 

Since consistency relation requires the stress point must remain 
on newly developed surface during isotropic hardening, we have 


dF 


8F 

3a ij 




+ 


9F • 



P 

ij 



(2.34) 


h can generally be expressed as a function of irrecoverable strain 
rates. 
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h . e p 


ae ?J « 


and substituting this in equation (2.34) yields 


dF - |L_; + . |£ J0-& = 0 

8<, iJ lJ 9 S Pj 1J 8h 9,Pj 1J 


(2.35) 


Equation (2.35) together with the normality condition (2.32) yields an 
expression for X 


X = 


L° 


ij CT ij 


where 


1 9F 


J ij 


H da. . 
1J 


(2.36) 


(2.37) 


H.-f + JL. (2.38) 

3 "mn 

Note the relationship given by equation (2.36) does not exist for 
ideal plasticity since H vanishes for zero hardening. This can be 
avoided by reformulating the above expression in terms of strain 
increments. In the absence of thermal strains we have 


where 



9F 




(2.39) 


(2.40) 


H 1 = 


6 or, 


iL. D e 
y klmn 


9F 


ki 


3 or, 


mn 


9F + 9F 9h 9F 
de kl ah 3cr kl 3or kl 


(2.41) 
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and it is obvious that for ideal plasticity H* does not vanish. In 
the above equations, an d are functions of the current state of 
stress. 

Having found a relation for X, the plastic strain rate can now be 
found utilizing the normality condition (2.32), and the elastic strain 
increment is determined using Hooke's law. 



— r P n 

- c ijkl CT kl 


(2.42a) 


5 ij = c ijkl ^kl 


(2.42b) 


where 


C p 

ijkl 



3F 


3 a. 


kl 


(2.43) 


The total strain increment is a summation of the plastic, elastic and 
thermal strain: 


e . . 
10 


pet 
8 ij + 8 ij + c ij 


where 


(2.44) 


ij 


8 ij * T 


(2.45) 


and therefore the mechanical strain can be defined as 


e ij = (e ij " ° 8 ij T) = 8 ij + 8 ij 


(2.46) 


or using equation (2.42) 


s m 

S ij 


C ep rt 

c ijkl ff kl 


(2.47) 


where 


C® p — p® + pP 

u ijkl ~ 4jkl + °ijkl 


(2.48) 
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Upon inversion of the above equation, a relation for the stress 
increment can be found in terms of the mechanical strain increment 


ep m 

°ij = D ijkl e kl 


(2.49) 


Where D e . p . is the elastop lastic constitutive tensor. More 

1. J X 

specifically for the isotropic, strain-hardening Von Mises material 


D?P 


3p S^4 S 


ij ~kl 




(2.50) 


2.5.2 Incremental Inelastic BEM Formulation 

The governing differential equation of equilibrium .is given in 
equation (2.11). Expressing it in incremental form: 


(X+n) u. .. + n u 

J l 


J i. jj 


+ f„- = 04 


lj .3 


(2.51) 


Boundary conditions, similar in form to equation (2.2), are 
imposed on an incremental scale and incremental body forces are 
applied in a manner analogous to section 2.4. 

The initial stress rate < 1 ?^ i s defined for the general case 
of thermoplasticity as a summation of the thermal stress rate and 
a plastic stress rate resulting from the nonlinearities present 
in the plastic domain. For the purpose of the thermoplastic BEM 
algorithm these can be defined as: 


0 t p 
®ij ~ c ij + c ij 


(2.52a) 


’ij 


(3X+2fi) oT 


(2.52b) 


a ij “ 


_e m ~ep m 
D ijkl 6 kl " D ijkl 8 kl 


(2.52c) 
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(2.52d) 


e kl = 8 kl “ 8 kl ° T 

where, »■ ^ijkl are e ^- as ^^ c an< * elastoplastic constitutive 

tensors respectively, E is the modulus of elasticity, a is the 

coefficient of thermal expansion, T is the incremental temperature 

change, and i S the mechanical strain. 

The axisymmetric, two-, and three-dimensional boundary integral 

formulation derived in section 2.2 and rewritten in incremental form 

are all applicable to inelastic BEM analysis. It is understood that 

the field variables are incremental quantities, and the nonlinearities 

are incorporated into the analysis through the initial stress rate c?. 

•J 

as defined in equation (2.52). 


2.5.3 Iterative Solution Algorithm for Thermoplasticity 

The algorithm described here provides the solution for an 
incremental, assembled system, analogous to the system defined by 

equation (2.27). The solution requires complete knowledge of the 

# 

initial stress distribution o° within the yielded region that is 
induced by the imposition of the boundary loading, body forces, and 
thermal loads. Unfortunately the nonlinear part of the initial 
stress is not known a priori for a particular load increment and 
therefore an iterative process must be employed within each loading 
stage. 

An important feature incorporated in the iterative algorithm of 
the present work is an iteration acceleration scheme (Raveendra, 1984) 
which utilizes plastic stresses generated by the past history. In 
this procedure, the path followed by the previous load increment is 
used to extrapolate the plastic stresses at the beginning of the 
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current increment before the iterative operations. This results in 
substantial reduction in computer time. This procedure is graphically 
illustrated in Figure 2.2 for a simple tension problem with a variable 
hardening parameter. The initial loading from point A generates 
plastic stresses which are distributed during the iterations to arrive 
at the solution at B . Upon loading from state B the extrapolated 
path BC is followed. At C the resulting corrective stresses are 
distributed to reach state D . 

The incremental algorithm is described below. Note, this 
algorithm requires that a stress equation be written at each cell 
node. 

(a) Obtain the elastic solution for an arbitrary increment of 

• • 
boundary loading y (and thermal loading T) from 

x = [A b ] -1 [B b y + cV) 

and (2.53) 

• • • • 
a = k a x + B°y + C 0 ^ 

(b) Scale the elastic solution such that the highly stressed node is 
at yield. In the case of nonproportional loading this onset of 
yielding cannot be reached by scaling. Instead, incremental 
loads must be applied until yielding is reached. 

• • 

(c) Impose a small load increment y (and T ), (usually less than 
five percent of the yield load) plus the estimated value of the 
plastic stress rates accumulated from the previous load step and 
evaluate 
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x = tA b ] _1 [B b y + C b 0°] 


(2.54) 


and 


= A°x + B^y 


C*0° 


where a° is a combination of the estimated plastic stress 
rates and the incremental thermal stress rates. If no prior 
plastic history exists, the value of the estimated 
plastic stress rates are zero. 


(d^ Accumulate all incremental quantities of stress, traction and 
displacement rates and use equation (2.52d) to calculate the 
mechanical strain of this increment. . 


(e) Evaluate the current constitutive matrix using the new stress 
history and calculate the current plastic stress rates via 



(D 


ijkl 


“ D ijkl^ 


*m 

e kl 


(2.55) 


Accumulate the plastic stress history of this load increment to 
be used as the estimated rate for the next load increment in step 
(c). 


(f) If the current increment of plastic stress rates, computed in 
equation (2.55), is greater than a prescribed tolerance (normally 
0.005 times the yield stress) at any node, then calculate the 
incremental quantities due to these rates using 

x = [Ah^tchoP] and a = k°x + cV* (2.56) 


39 



and return to step (d) for the next iteration. If the value is 
less than the prescribed tolerance go to step (g). Note the 
boundary and thermal loading is zero for this calculation and the 
mechanical strain rate will therefore be equal to the total 
strain rate. If the number of iterations is greater than a 
specified limit (usually 50) the system is assumed to have 
reached the state of failure. 

(g) Return to step (c) and apply the next load increment and the 
accumulated plastic stress rates from this load step. (If the 
size of the load increment changes the estimated plastic stress 
rates should be scaled proportionally.) 

Any residual plastic stress at the end of the iteration is carried 
forward and applied to the system with the next load increment. 

2.5.4 Variable Stiffness Plasticity Approach 

The above nonlinear formulations include initial stress rates in 
the governing equations which are not known a priori and therefore are 
solved by using iterative procedures. A non-iterative, direct 
solution procedure is made feasible by reducing the number of unknowns 
in the governing equations by utilizing certain features of the 
incremental theory of plasticity expressed by equations (2.32), (2.3 6 ) 
and (2.3 9). This non-iterative or variable stiffness approach was 
introduced by Raveendra (1984) in a single-region, two-dimensional 
context. For the first time, this method is implemented in an 
axisymmetric , two— and three-dimensional, multiregion computer 
program. 
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The initial stress rates cr?^ appearing in the incremental form 
of integral equations (2.20) to (2.22) can be expressed in the context 
of an elastoplastic deformation as: 


ff ij " K ij * 


(2.57) 


where K ±J = D? jkl 


8F 


3cr, 


kl 


Substituting equations (2.36) and (2.57) in equation (2.20) and (2.22) 
we obtain (in the absence of body forces): 

C ij u i(?) = | [G i j(x,^)t i (x) - F i j(x,?)u i (x)]dS(x) 


+ J B ± j(x,5)K ± (xU(x)dV(x) 


(2.58) 


and 


i ( 5> " ^ k <?> J [aJ Jk (x,{)t i (x) - F?j k (x,^) u i(x)]dS(x) 

S 


ijk' 


+ L j k (?) J Bj pJk (x,§)K ip (x)X(x)dV(x) 
V 


(2.59) 


Equations (2.58) and (2.59) can be solved simultaneously to evaluate 
the unknown values of displacements, traction rates and the scalar 
variable \ . 

Although equation (2.59) can be applied to any elastoplastic 
strain hardening problem, becomes indeterminate for the case of 
ideal plasticity (zero strain hardening), therefore, it cannot be used 
for ideal plasticity. However, this minor problem can be circumvented 
either by specifying a small amount of strain hardening or more 
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appropriately replacing equation (2.59) by the strain rate equations 
(2.39) and (2.21). Thus using equations (2.39) and (2.57) in (2.21) 
we obtain 

X<$) = L jk (5) I [G ijk (x ' ?)t i (x) " F ijk (x ' d)u i (x)]dS(x) 

S 

+ L j k U) J B® pjk (x.?)K ip (x)X(x)dV(x) (2.60) 

V 

Equation (2.60) can now be applied to any problems of ‘elastoplastdcity 

(both strain hardening and ideal plasticity). 

Assembly of System Equations for the Variable Stiffness Method - 

Although equations (2.36), (2.3 9) and (2.57) were substituted into the 

analytic form of the integral equations for the presentation, during 

actual implementation, these equations are substituted, on a nodal 
* 

basis, into the incremental form of the assembled equation system 

(2.27). Therefore, using the incremental form of equation (2.27) as a 

starting point, equations (2.36) and (2.57) are employed to arrive at 

the following system: 

* # • 

A b x = B b y + C^X 

and (2.61) 

X = LA®x + LB®y + LC®KX 

. • 

where y are the known incremental boundary conditions and x are the 

unknown. Matrices A, B, C and A®, B®, C® are constant for a given 
problem and matrices K, L are dependent upon state variables, and 
are assumed to be constant during each small load step. 
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It should be noted that the second equation is written only for 
the cell nodes that are expected to yield during the current load 
step. 

The above equations can be rewritten as: 

A b x = b b + C Xb X (2.62a) 

and 

X - A** + b X + C X X (2.62b) 

or upon rearranging the equation (2.62b) we have 

HX = A*x + b X (2.63) 

• • 

where b b = B b y , C xb = C 1 ^ , A x = LA 0 , 

b x = LB 0 y , C X = LC?Z , H = I - C x , and 
I is the identity matrix. 

Equation (2.63) can then be recast as 

X = A°i + b° (2.64) 

where 

A° = H" 1 A X 
b° = H" 1 ^ 

Substituting the above equation into equation (2.62a) results in the 
final system equations: 

A r x = b r (2.65) 

where 
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A r = A b - C Xb A° 


b r = b b + C Xb b° 

• 

The above equation can be solved to evaluate the unknown vector x at 
boundary nodes for every increment of loading. The present 
formulation is similar to the variable stiffness approach used in the 
finite element method since the system matrix on the boundary as well 
as the right hand side vector are modified for each increment of 
loading. For a multi-region system the inversion of H can be carried 
out for each region separately. 

Solution Process for Variable Stiffness Approach - As previously 
mentioned, the solution process does not involve any iterative 
procedure, instead the substantial part of the solution effort is 
spent on assembly of the system equations for each load step. These 
operations can be described as follows: 

(a) Impose an arbitrary boundary loading and solve the elastic 
problem in the usual manner. 

(b) Scale the elastic solution such that the highest stressed node is 
at yield. 

(c) Apply a small load increment (usually < five percent of the yield 
load) and compute K and L matrices using the past stress 
history. 

(d) Form the system equation (2.65) and solve for x. 

(e) Evaluate the initial stress rates <y° using equations (2.57) and 
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(2.64) : 


cr° = K A = K[A°x + b°] 

(f) Evaluate interior quantities displacement and stress rates using 
the incremental form of equations (2.24) and (2.25). 

(g) Return to step (c) if the strains are less than a specified 
norm. Otherwise, failure is assumed to have occurred. 

It is of importance to note that the matrices K, L do not exist in 
the elastic region, therefore, the corresponding equations involving 
these matrices are formed and a° is determined, only for the nodes 
that are at yield. Any small deviation from the yield surface can be 
corrected by applying the stress rate difference (i.e. the initial 
stress rate) during the next load step. 

2.6 ACCURACY AND CONVERGENCE 

In order to test the accuracy and convergence of the 
elastoplastic formulations, simple problems with known solutions, such 
as cubes, spheres, and cylinders, were analyzed. Axisymmetric, two- 
and three-dimensional inelastic formulations were coupled with both 
the iterative and variable stiffness solution algorithms. All 
possible combinations were tested and in all cases excellent agreement 
with the analytical solutions was obtained. A few of these examples 
are presented here. 

2.6.1 Two-dimensional Analysis of a Thick Cylinder 

In order to verify convergence and accuracy of the two- 
dimensional. iterative plasticity algorithm, an analysis of a thick 
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cylinder subjected to increasing internal pressure is carried out 
under the plane stress condition. The discretization of the cylinder 
(1:2 ratio), shown in figure 2.3, has ten elements and four volume 
cells. In figure 2.4, the hoop strain vs. pressure for the inner and 
outer surfaces are compared to the analytic solution (Hill, 1950). 
Excellent agreement is obtained. Other numerical results, to within 
1% error of the analytical solutions were obtained for the radial and 
hoop stress through the cylinder at different load pressures, and the 
displacement versus internal pressure response. In figure 2.5, the 
axial strain vs. pressure is compared with the theoretical solution, 
and slight deviation is observed. A possible explanation for this 
error is that the theoretical solution assumes the Tresca yield 
criterion, where the present analysis is based on the Von Mises 
criterion. The computational time using the two-dimensional analysis 
was found to be approximately 5% of the time required for the three- 
dimensional BEM analysis. 

2.6.2 Three-dimensional Analysis of a Thick Cylinder 

A three-dimensional elastoplastic analysis of a thick cylinder 
(in plane strain) is carried out using the variable stiffness 
algorithm. The discretization of the cylinder (1:2 ratio), shown in 
figure 2.6, has eighteen boundary elements and four (twenty-noded) 
isoparametric cells. The material properties are (given in consistent 
units) : 

E = 2600 
v =0.3 

<* 0 = 600 (Von Mises yield criterion assumed) 
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In figure 2.7, the load-displacement response (assuming ideal 
plasticity) is compared to the theoretical result (Hill, 1950). In 
figure 2.8, the circumferential stress distribution through the thick 
cylinder is presented for two different load levels. Good agreement 
with the analytical solution is observed in both figures. 

In a second analysis, the material is assumed to exhibit variable 
plastic strain hardening according to: 


Stress vs. 

Plastic 

600 

0.0 

640 

0.1 

660 

0.2 

660 

10.0 


In figure 2.9, the load-displacement response at the inner and 
outer surface is presented. Results obtained from the iterative 
procedure (with 5% load increments) are compared with two sets of 
results from the variable stiffness method. One set is obtained using 
2% load increments, and the other 5% increments. Results are 
generally within 1% of one another. 

2.6.3 Axisymmetric Analysis of a Hollow Sphere 

A hollow sphere (1:2 ratio) subjected to internal pressure is 
used to demonstrate the axisymmetric elastoplastic formulation. Nine 
boundary elements and ten (eight noded) cells are used in the 
discretization of the sphere as shown in figure 2.10. The mesh 
utilizes spherical symmetry and contains points on the origin. The 
incline face is subjected to a roller boundary condition. The load- 
displacement response is shown in figure 2.11 and the hoop stress 
through the sphere at different load levels is given in figure 2.12. 
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The BEM solution obtained using the iterative method is in good 
agreement with the analytical result (Hill, 1950). 

2.7 CONCLUDING REMARKS 

The conventional boundary element formulations for elastic and 
inelastic thermal stress analysis were presented. The inelastic 
axisymmetric implementation presented in this chapter, with its use of 
quadratic elements and multi-region facility, is the most advanced 
implementation of its kind. Two inelastic algorithms, an iterative 
and a variable stiffness type approach were employed, the latter for 
the first time in axisymmetric and three-dimensional analysis. By 
comparing results of test problems to their analytical solutions, the 
accuracy of both methods was demonstrated. In Chapter 6, the methods 
will be used to analyze larger problems of practical interest. 

The body force and nonlinear effects in the present formulation 
were incorporated -in the system through volume integrals or extra 
surface integrals. In the next two chapters, new approaches for the 
treatment of these effects will be presented. 
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table 1 

STRESS SmiES FOR INITIAL STRESS EXPANSION 
IN TWO-DIMENSIONAL PLANE STRAIN 
(PLAIN STRESS) ANALYSIS 


Coefficient to 

Stress State be determined 
corresponds to 


where 

E 


Nodal Values of Assumed Stress State 

XX 


•Sr 

“x 
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0 

0 

X 

? 

1 

vH 

-v(l+v)y 

0 

E 

0 

-v(l+v)x 

(l-v 2 )y 

0 

0 

E 

(l+v)y 

(l+V )x 


is the modulus of elasticity, 
is the Poisson's ratio, and 


x and y are the nodal coordinates. 

All stresses and tractions are zero for all stress states. 

The stress states for two-dimensional plane strain analysis, given in the table 
above, can be applied to the plane stress case, if the modified material parameters, 
defined belcw, are used. 

- ^ E(l+2*) 7 - — 

(l+v) 2 1+V 






























TABLE 2 


STRESS STATES FOR INITIAL STRESS EXPANSION TECHNIQUE 
IN THREE-DIMENSIONAL ANALYSIS 


Stress State 

Coefficient to 
be determined 
corresponds to 

Nodal values of Assumed Stress State 

CT xx 

°?y 


„0 

xy 

4z 

■Si 



°z 

1 

_o 

°xx 

E 

0 

0 

0 

0 

0 

X 

-vy 

-\>Z 

2 

°YY 

0 

E 

0 

0 

0 

0 

-vx 

y 

-\>Z 

3 
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0 

E 

0 

0 

0 

-VX 

-vy 

Z 
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0 

0 

0 

E 

0 

0 

(l+v)y 

(l+v)x 

0 

5 

°Zz 

0 

0 

0 

0 

E 

0 

(l+v)z 

0 

(l+v)x 

6 


0 

0 

0 

0 

0 

E 

0 

(1+ v)z 

( 1+ M)y 


where 

E is the modulus of elasticity, 

v is the Poisson's ratio, and 

x.y, and z are nodal coordinates. 

All stresses and tractions are zero for all stress states. 



TABLE 3 


STRESS STATES FOR INITIAL STRESS EXPANSION TECHNIQUE 
IN AXISYHHETRIC ANALYSIS 


Stress State 
number 

Coefficient to 
be determined 
corresponds to 

Nodal values of assumed stress states 

*rr 

a zz 


°rz 

o 

°rr 

o 

°zz 

0 

0 

°rz 

u r 

U 2 

fc r 
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1 
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rr 
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r r 
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c rr n r 

Vz 

3(l-v)C r 

(l-v)C * 

3(l-v)C 1 

2 

o 

a rr 
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E 

0 

E 

0 

(l-v)r 

-2vz 
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zz 
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0 

0 
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— vr 
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■ 

4 

o 

a rz 

0 

m 

■ 

■ 

H 

0 

0 

■ 

0 

■ 

■ 

■ 


where E modulus of elasticity 

fi shear modulus 

v Poisson's ratio 


r and z are nodal coordinates 

n and n are normals of the boundary 
r z 

C arbitrary parameter with dimensions of length. added to insure dimensional homogeneity. 

The value can be equated to one for simplicity. 







Figure 2.1 

Isotropic Hardening Behavior 
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Figure 2.2 

Accelerated Iterative Scheme 
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Figure 2.4 

Hoop Strain in a Thick Cylinder under Internal Pressure (Plane Stress) 
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Figure 2.5 

Axial Strain at the Outer Surface of Thick Cylinder under 
Internal Pressure (Plane Stress) 
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Figure 2.7 

Radial Displacement of the Outer Surface of a Thick Cylinder under 
Internal Pressure (Plane Strain) 
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Figure 2.8 

Circumferential Stress Distribution through a Thick Cylinder under 
Internal Pressure (Plane Strain) 
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Figure 2.9 

Load versus Radial Displacement of the Inner and Outer Surfaces of a 
3-D Thick Cylinder with Hardening (Plane Strain) 




Figure 2.10 

Hollow Sphere (1:2 Ratio) 
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Figure 2.11 

Radial Displacement in a Hollow Sphere subjected to Internal Pressure 




Figure 2.12 

Hoop Stress in a Hollow Sphere subjected to Internal Pressure 
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CHAPTER THREE 


ELASTIC BOUNDARY ELEMENT FORMULATIONS 
BASED ON PARTICULAR INTEGRALS 


3.1 INTRODUCTION 

In this chapter, a novel boundary element approach for the 
treatment of body forces is presented. The method is unique since no 
volume integration or additional surface integration is required as in 
earlier formulations. The method is based on the well known concept 
of developing the solution of an inhomogeneous differential equation 
by means of a complementary function and particular integral. 

In the first section, the particular integral formulation for a 
general inhomogeneous differential equation is developed and the 
boundary element solution procedure is presented. The method is then 
applied to gravitational and centrifugal body force analyses. 
Finally, particular integrals are developed for thermal body forces 
according to the uncoupled theory of thermoelasticity. In this theory 
the heat conduction problem is independent of stress, and therefore, 
will not be considered here. Instead, the temperature distribution is 
assumed to be a known (or previously solved) function. 

The axi symmetric, two-, and three-dimensional formulations are 
implemented in a general purpose, multi-region system, and examples 
are presented to demonstrate the accuracy of the method. 
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3.2 BOUNDARY ELEMENT FORMULATION USING PARTICULAR INTEGRALS 

3.2.1 General Theory 

A solution satisfying a linear, inhomogeneous, differential 
equation and boundary conditions can be found using the method of 
particular integrals if the complete solution of the corresponding 
homogeneous equation is known, provided, of course, a particular 
solution can be found. The procedure is described below: 

A linear inhomogeneous differential equation for a vector 
variable u.^ can be expressed in operator notation as 

m»i> + fi - 0 

where 

L(u^) is a linear, self-adjoint, ' homogeneous, differential tensor- 
operator, and 

f^ is the inhomogeneous function. 

The solution of the above equation can be represented as the sum 
of a complementary function u® satisfying the homogeneous equation 

L(u?) = 0 (3.2) 

and a particular integral u|? satisfying the inhomogeneous equation 

L(u?) = f i (3 - 3) 


% 

(3.1) 


The total solution u^ is expressed as 


U ± = u i + ul 


(3.4) 
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In the theory of linear, inhomogeneous, differential equations, 
it is understood that the particular integral is not unique, and any 
expression satisfying equation (3.3) is a particular integral, 
regardless of boundary conditions or how it was obtained. 

Particu l ar Integral - The particular integral is classically 
found via the method of undetermined coefficients, the method of 
variation of parameters, or obtained by inspection of the 
inhomogeneous differential equation. When these techniques fail to 
produce a particular integral, a general procedure can be used. Here, 
the singular solution satisfying the homogeneous differential equation 
is multiplied by the inhomogeneous quantity, and the product is 
integrated over an infinite domain. The final result is a particular 
integral of the inhomogeneous equation. This can be expressed 
mathematically as 

00 

«?U> = f G..(x,£) f,(x) dV(x) (3.5) 

x ico AJ J 

where G. .(x,£) is the fundamental singular solution of the homogeneous 
differential equation. The use of a polar coordinate system will 
simplify the integration. 

Once a particular integral is found, it is added to the 
complementary function to form the total solution. The parameters of 
the complementary function are adjusted to insure that the total 
solution satisfies the boundary condition, and hence, produces a 
unique solution to the boundary value problem. 
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з. 2.2 BEM Solution Based on Particular Integrals 

The use of particular integrals in BEM was tentatively discussed 
by Watson (1979) and Banerjee and Butterfield (1981). In 1986, Ahmad 
and Banerjee (1986) successfully employed the idea in a two- 
dimensional eigenvalue analysis. Gravitational and centrifugal 
formulations have been presented in axisymmetry (Henry, Pape, 
Banerjee, 1987) and .two-dimensional (Pape, Banerjee, 1987) analysis. 

Equation (3.1) can be defined for an elastostatic stress analysis 
with gravitational, centrifugal or thermal body force as 

L(u i ) + f ± - 0 (3 ‘ 6) 

where 

L(u i ) = u j,ij + 

и, represents the displacement vector, 
i 

f^ = -pge^ for gravity loading in the z direction, 

f^ = pco^(x^e 2 + X 2 e 2 ^ f° r centrifugal loading about the z axis, and 

f^ = -pT,^ for thermoelastic loading 

When more than one type of body force is present, the net body force 

is the summation of the individual body forces. 

Complementary Function - The boundary integral equation 
satisfying the homogeneous part of the differential equation is the 
complementary function in this procedure. The complementary function 
for displacement at point is expressed as 

C ij ( 5 ) u i ( 5> = I [Gij(x,5> tj(x) - F^tx,?) uj(x)] dS(x) 

* (3.7) 

where the u^ and t^ are the complementary functions for displacement 
and traction, respectively. The total solution for displacement u^ 
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and traction ^ are 
* u° + u£ 

\ - tj + tj (3.8) 

where u p and t? are the particular integrals for displacement and 
traction, respectively. 

Explicit expressions for u p and t^ will be derived for 
gravitational, centrifugal and thermal body forces in subsequent 
sections, but before proceeding to the formulation of these particular 
integrals, a brief overview of the particular integral based BEM 
procedure is presented. 

Method of Solution - The boundary integral equation (3.7) is 
discretized and integrated for a system of boundary nodes in the 
manner described in Chapter 2. The resulting equations can be 
expressed in matrix form as 

G t° - F u c =0 (3.9) 

By introducing equation (3.8) into equation (3.9), the 
complementary function can be eliminated; 

Gt - F u = G t p - Fu p (3.10) 

The particular integral terms on the right hand side of this equation 
are functions of known body forces. 

In a multi-region system, a set of (complementary) equations, 
similar to equation (3.9) are generated independently for each region. 
Likewise, the particular integrals of each region are derived 
independently, leading to a set of equations for each region similar 
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in form to equation (3.10). Interface conditions, expressing the 
interaction of real quantities between region, are applied, and after 
assembling the unknown boundary quantities and corresponding 
coefficients on the left hand side, and the known boundary conditions 
on the right, the final system can be written as 

A b x = b b + b p (3.11) 

where A b is a fully populated matrix, vector x represents the unknown 
boundary conditions, vector b b is the contribution of the known 
boundary conditions and vector b p is the contribution of the 
particular integral. This equation system can be solved for the 
unknown vector x by standard numerical techniques. 

A Note on Multi-region Programming - The most efficient solution 
procedure for a single region BEM system utilizes equation (3.9) 
instead of equation (3.10) since equation (3.10) requires additional 
matrix multiplication. In this procedure the particular integrals u^ 
and t p are calculated and the complementary functions corresponding to 
the known boundary conditions are determined using 

u i = U i (x) " u i and b i = T i (x) " fc i 

where UNfx) and T^x) represent the imposed displacement and traction 
boundary conditions, respectively. The solution to equation (3.9) 
yields the unknown complementary functions, and when added to the 
corresponding particular integral the correct solution is obtained. 

This procedure is not as straight forward for multi-region 
implementations, since interface conditions relate real quantities of 
displacement and tractions, not the complementary quantities. In a 
large, multi-region system with a complex assembly for local 
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definition of boundary conditions, sliding between interface elements, 
etc., it may be beneficial to use equation (3.10) in spite of the 
extra matrix multiplications rather than equation (3.9), since the 
latter would require extensive modifications of interface conditions. 

3.3 PARTICULAR INTEGRALS FOR GRAVITATIONAL AND CENTRIFUGAL 

BODY FORCES 

In the previous section, a boundary element formulation utilizing 
particular integrals was described. Application of this procedure 
requires the determination of the particular solution for displacement 
and traction at all boundary nodes. Particular integrals for 
gravitational and centrifugal body force loadings, applicable to two- 
dimensional, three-dimensional and axisymmetry boundary element 
analysis, are presented in this section. 

The two- and three-dimensional particular integrals are most 
readily obtained (Sokolnikof f , 1956) by the method of undetermined 
coefficients. Using tensor transformations, the axisymmetric 
particular integrals can be obtained from the three-dimensional form 
(Henry, Pape, and Banerjee, 1987). For pure axisymmetry, i.e., 
axisymmetric bodies under axisymmetric loads these transformations 
simplify to the following: 

Coordinates: Xj = r 

x 2 = o (3.12) 

x 3 “ 2 

Displacement: u r = Uj 

u 9 = 0 (3.13) 

u z - u 3 
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Stress: 


rr 

= <r n 

a 00 

= °22 

a zz 

= c 33 

a rz 

= a 13 

a r0 

II 

N 


and tractions are found using the 


(3.14) 


Cauchy traction relation 


fc i = »ij n j i ’ i = r ' Z i . * 

The application of the above transformation are straight forward and 
therefore only the two-, and three-dimensional particular integrals 
will be defined. 

The following two— dimensional , particular integrals assume the 
plane strain conditions. 

The two-dimensional particular integrals for the plane stress 
condition are obtained from the plane strain conditions replacing the 
material constants v and X by a modified constants v and X 


v = v/(i+v) 


(3.16) 


7 = (3.17) 

K X+2(i 

The shear modulus n remains unchanged. 

3.3.1 Particular Integrals for Gravitational Body Forces 
T wo-dimension (plain strain) - The body force in a two- 
dimensional solid of density p under gravitational acceleration g 
directed along the -X 2 axis is expressed as 
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f l = 0 


(3.18) 


f 2 = “PS 

The particular integral for displacement under this body force is 


u 


P _ Pg 
1 4p(X+p) 


Xxix 2 


U P ~ pg . 

2 8 t i(X+ ( t) 


[(X + 2(i)x 2 + Xx^l 


(3.19) 


and the corresponding particular integral for stress is 


a 22 = "P6 x 2 


°11 “ *12 “ 0 


and the particular integrals for tractions is 


(3.20) 




x 2 


n 2 


(3.21) 


Three-dimension - The body force in a three-dimensional solid of 
density p under gravitational acceleration g directed along the -x^ 
axis is expressed as 



(3.22) 


The particular integral for displacement under this body force is 


,P = ^£8 

1 2p(3X+2(i) 


x l x 3 
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p Xpg 

2 " 2|i(3X+2p) X 2 X 3 


(3.23) 


„p _ -^ + P)Pg v 2 _ *Pg _ (x 2 + 2) 
u 3 “ 2p(3X+2p) x 3 4 |a(3X+2|i) ^ X 1 x 2 


and the corresponding particular integral for stress is 


- "Pg x 3 


°22 ~ ^12 + ct 13 


and the particular integral for traction is 


= t£ = 0 


— pg x 3 n 3 


(3.24) 


(3.25) 


3.3.2. Particular Integral for Centrifugal Body Forces 

Two-dimension (plain strain) - The centrifugal body force in a 
two-dimensional solid of density p rotating at a speed w about an axis 
parallel to the x^ axis and centered at a reference coordinate is 


expressed in indicial notation as 


> ** , r 
i = p«o y± 


(3.26) 


where 


y i = x i " «i 

The particular integral for displacement under this body force is 


U i " 8(X+2p) (y k y k> x i 


(3.27) 


where summation is implied over k. 
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The corresponding particular integral for stress and traction is, 
respectively 


•ij- 


4<!w2m> t<2X+|,) Vk 6 ij * 2 K 


(3.28) 


fc i " 40Sl [(2X+tl) y k y k n i + ^ y i y k n k ] 


(3.29) 


where 

is the boundary normal at x^, and 
is the Kronecker delta. 

Three-dimension - The centrifugal body force in a three-dimension 
solid of density p rotating at a speed w about an axis parallel to the 
Xj axis and centered at a reference coordinate is expressed as 


f 

f 

f 


1 

2 

3 


2 

= pt» Yj 
2 

= pu y 2 
= o 


where 


(3.30) 


y i = x rh 

Note, ^ = 0 in axisymmetry. 

The particular integral for displacement 


under this body force is 


,p _ 

-p(j)^ 1 

r 5X+4p 

<y i^2> + y l 1 

i - 

8(X+2|i) 1 

L 4(X+p) 

p _ 

-pu> 2 I 

r 5X+4p 

(y l* y !> + 55f *3 3 

2 “ 

8(X+2|i) 1 

L 4(X+p) 
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(3.31) 



2+3 V+ 2 V 2 -1+5 V +2^ 2 -( 2 +v) (l-2v) 

C 2 = 16(l-v) °4 = 8(l-v) °6 " 8(l-v) 

and the particular integral for traction is obtained using the Cauchy 
traction relation 

t jL = nj i» j = 1,2,3 (3.33) 

3.4 PARTICULAR INTEGRALS FOR THERMAL ANALYSIS 

The gravitational and centrifugal particular integrals, presented 
in the previous section, were relatively straight forward since their 
body forces had specific distributions that were represented by 
elementary functions. In thermal stress analysis this is not the 
case, since thermal body force can exhibit general temperature 
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distributions. So consequently, thermal particular integrals are more 
difficult to obtain. 

In this section, the aim is to develop a particular integral 
applicable to any temperature distribution. This distribution is 
expressed by nodal temperatures specified at random point through the 
domain of the body. 

The present formulation starts with the assumption that the 
particular integral for displacement can be expressed as a gradient of 
a thermoelastic displacement potential h(x) in accordance with the 
development of linear quasi-static thermoelastic theory (Nowacki, 
1962) : 


u^(x) = k h.^x) 


(3.34) 


where 


q( 3X+2n) 
U+2n) 


(3.35) 


Substituting equation (3.34) into equation (3.6) and simplifying 
yields 


h,jj(x) = T(x) (3.36) 

Hence, an implicit relation between u£ and T is derived with the 
use of a potential function h(x). The irrotational form of the 
particular integral in equation (3.34) does not imply a loss of 
generality in the total solution since the balance of the temperature 
effect is incorporated into the complementary function. 
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The particular integral of equation (3.34) is of little interest 
unless a suitable function h(x) can be found to accurately represent a 
general temperature distribution in equation (3.36). This is achieved 
using a device, first introduced in the context of finite elements, 
called the global shape function. 

3.4.1 Particular Integrals for Two- and Three-dimensional 
Analysis 

Assuming the function h(x) can be represented by an infinite 

V ft 

series, an expression relating h(x) to a set of fictitious scalar 
densities 0(? n ) via a global shape function C(x,£ n ) can be written as 


h(x) - J C(x,5 n )0<5 n ) 
n=l 


(3.37) 


where C(x,£ n ) = a' suitable function of spatial coordinates x and 
^n* 

Several functions were considered for this purpose, however, the 
best results were obtained with the following expression: 


C(x,? n ) , A 2 tp 2 - b n p"] 


(3.38) 


where 


A is a characteristic length, 
o 


p is the euclidean distance between the field point x and the 


source point § , and 


b R is a suitably chosen constant to be discussed later. For 

the present discussion, assume b n ” 

All distances are nondimensionalized by a characteristic length 
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Using equations (3.34) and (3.37), the particular integral for 
displacement is found: 

00 

uj(x) = J D i (x,5 n )0(§ n ) (3.39) 

n=l 

i * 1,2 for two-dimensions 
i = 1,2,3 for three-dimensions 

where 

D j.(x.5 n > - HC.jlx.?,,) = kA 0 [2 - 3b n plyj 
y i * '*1 - «„>i> 

and 

d = 3 for three-dimensional or d = 2 for two-dimensional ' (plane 
strain) analysis. 

Applying the Laplaeian operator on equation (3.37), an expression 
for the temperature distribution is derived 


00 

T(x) = J K(x,f )0<£ n ) (3.40) 

n=l n n 

where 

K(x,? n ) = C. ii (x,§ n ) 

= [2d - 3(l-*i> b p ] 
n y 


A particular integral for strain can be found upon substitution 

of equation (3.39) into the strain-displacement relation. 

00 

e kl ( *> • I E kl <*.? n )0(5 n > <3 - 41 > 

n=i 

where 


E k l(x,5 n ) = kC, kl (x,? n ) 


y k y l 

= k [28 kl - 3b n (8 kl p + -f±)l 
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This result introduced into the stress-strain law of 
thermoelasticity produces the particular integral for stress: 

09 

- X h 3 (*.?„>»<{„> < 3 -«> 

where 

S ij ^’^n 5 “ D ijkl E kl (x ’^n^ “ 8 i j and 

D ijkl = x 5 ij 6 kl + 2v - 8 ik 8 jl ' p = o(3X+2fl) 

Finally a particular integral for traction is derived by 
multiplying the above equation by appropriate normals 

CO 

tj(x> = \ H i (x,5 n )0(? n ) (3.43) 

n=l 

where 

H i(x,5 n ) = Sij(x,5 n ) n j(x) , and 

n.(x) = unit normal at x in the direction. 

Once again, the plane stress formulation can be obtained from the 
plane strain case by substituting the modified material constants o 
and X into the plane strain equation in place of o and X, 
respectively, where 

_ o(3X+2fi) 

a “ 2(2X+|i) 

and X is defined in equation (3.17) 

3.4.2 Particular Integrals for Axisymmetric Analysis 

The axisymmetric particular integrals can be derived in one of 
the two ways described below: 
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1 . 


Direct Derivation - The particular integrals in equations 
(3.34) and (3.36) are expressed independent of coordinate 
systems using operator notation. Upon choosing a suitable, 
cylindrical form of a global shape function C(x,$), the 
axisymmetric particular integrals are then obtained by 
applying the axisymmetric form of the differential 
operators. Note, a suitable shape function must render a E 
matrix that is well-conditioned and will admit to an 
inversion (see equation 3.53). 

2. Indirect Derivation - The three-dimension particular 
integrals, given in equations (3.39) through (3.43), are 
rewritten in cylindrical coordinates. Upon application of 
appropriate tensor transformations and analytic integration 
in the angular (U) direction, a suitable set of axisymmetric 
particular integrals are obtained. 

The resulting particular integrals should exhibit axisymmetric 
behavior and satisfy the physical conditions at the origin (u£=0, 

e r Z =0» 8 rr“ 8 &(P* Choosing a global shape function for the direct 
derivation that will satisfy these conditions and still admit to an 
inversion is a difficult task. Nevertheless, the axisymmetric 
particular integrals derived from the three-dimensional functions 
using the indirect derivation should inherently satisfy these 
conditions. For this reason, the indirect derivation has been adopted 
in this dissertation. The resulting axisymmetric particular integrals 
are given below. 
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The particular integral for displacement* temperature and strain 
are, respectively: 


uP(x) - Jd^x,^) 0( V 

n=l 


( 3 . 45 ) 


i = r,z 


T(x) = ^ K(x,§ ) 0<5 n > 
n=l 


( 3 . 46 ) 


8 ?1<*> ®<*n> 


( 3 . 47 ) 


ij = rr,zz,rz,©0 


where D^x,^), K(x,5 n ), E ± j(x,$ n ) are subsequently defined. 

In the present formulation, strain is derived from displacement 
where differentiation (with respect to the field point x) is performed 
prior to the analytic integration (with respect to the source point 
? n ) in the angular direction. 

The particular integral for stress and traction are derived 
directly from strain in the usual manner employing the following 
relations 


°ij " D ijkl 8 kl 8 ij 8 T 


( 3 . 48 ) 


fc i " °ij n j 

The following notation is used in defining the axisymmetric 
particular integrals. 


r x = x r rad;1 - al coordinate of the field point 


z = v axial coordinate of the field point 
*x z 


r e = (F ) radial coordinate of the fictitious density node n 
q n 
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z £ = U z ) n axial coordinate of the fictitious density node n 

Note, the subscript n is dropped for simplicity. All coordinates 
are nondimensionalized by a characteristic length k Q . 

z = z x - 

R = [(r x + r^) 2 + z 2 ] 1/2 

_ a( l± vj 
K " (1-v) 

m = 4 v ?/ r2 

= o-l 

K = K(m) is the complete Elliptic Integral of the first kind 
E = E(m) is the complete Elliptic Integral of the second kind 
E and K are defined in Appendix III.D 
General Form (r x £ r £ and/or z x t z £ ) : 

°1 =K 

o 2 = (E-n^lO/m 

c 3 = [2E(2m-l)+m 1 K(2-3m)]/3m 2 

F 1 = 4c i /R 

F 2 = -4(2c 2 ~c 1 )/R 
Fj — 4(40j-4c 2 +c 2 )/R 
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Oj = 4RE 


0 2 = 4mR(cj-C2> 

If r =o and/or r.=0 and z v tz f , then the following substitutions should 

x s * 5 

be made for F and G„ 
a a 

F 1 = 2n/R 
F 2 " 0 

L 4 * 

p 3 = Fj/2 
G 1 = 2jtR 
g 2 = o 

Particular Integral for Temperature: 

K(x,£ n ) = 2n - 3b n Gj 

Particular Integral for Displacement: 

0 r (x,5 n ) - kA 0 [f 1 r x - | b n (r x G 1 -r ? G 2 )l 

V X >V - l F " I b n G l ! 

Particular Integral for Strain: 

E rr<*-V * * <r - ! b n < V-'xV^V'IV 1 

*„(*.«.) -k 'r-! b n <Gi+*fri>] 

E rz (x '*n ) = kz t- | b n (r x F r r 5 F 2 )] 
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E ee (x '^n ) = D r (x '?n )/(r x A o ) 

If r x =0» then 

*««<*• V * E rr <3 ‘-«n> 

Singular Form (r y = r g and z x - z g ) : 

The general form of the axisymmetrie particular integrals are 
singular when r x = and z x = z^. In this circumstance, the limiting 
form of these functions, given below, should be used. 

Particular Integral for Temperature: 

K(x,£ n ) = 2n-24 b n r x 

Particular Integral for Displacement: 

- k# o r * [ r - 8 Vx 1 

=0 

Particular Integral for Strain: 

E rr (x>£ n ) = k [3 10 b n r x ] 

E zz^ x '^n* = k ^3 ^ b n r x^ 

E rz (x 'V = 0 

E e9 (x,5 n ) - k t|i - 8 b n r x ] 
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3.4.3 Numerical Implementations 

The particular integral for displacement u£(x) and traction t£(x) 
must be evaluated at each boundary node before a solution to the 
governing equation (3.6) can be achieved. For the purpose of 
numerical evaluation, the infinite series representations of the 
particular integrals are truncated at a finite number of N terms. 
(The choice of N will be discussed later.) The particular integrals 
for displacement, traction, and temperature distribution can be 
rewritten for N number of terms as follows: 

N 

U i (x) = 5 Di(x,5 n )0<5 n ) (3.49) 

n=l 

N 

t£(x> = J H i (x,5 n )0(5 n ) (3.50) 

n=l 

N 

T(x) = ^ K(x,5 n )0($ n > (3.51) 

n=l 

The particular integrals are derived for each region, independent 
of the other regions. Hence a set of (the above) equations are 
written for each region where a temperature change is prescribed. 

The evaluation of the uP and t£ in the first two equations 
requires the determination of N fictitious nodal densities 0(5 n )» n = 
1 to N. For this reason N temperature equations (3.40) or (3.46) are 
written at each g n node. These equations are expressed in matrix 
form, independently for each region, as 

T = K 0 (3.52) 
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where K is a N x N matrix. Since the increment of temperature 
distribution is known, the fictitious nodal values 0(£ n ) can be 
determined: 

0 - K _1 T (3.53) 

A back substitution of the quantities 0(5 n ) in the particular 
integrals for displacement and traction allows for their evaluation at 
any point. Hence, u^ and t^ can be determined at all boundary nodes, 
and the boundary value problem can be solved as discussed in section 
3 .2.2. 

Disp l acements, Stresses, and Strains at Interior Points- The 
solution for the interior quantities consist of a complementary and a 
particular part: 

u i ■ u i + u i 

*lj * e ij + e iJ <3-54> 
c p 

CT ij = "ij + ff ij 

The particular integrals of the above quantities are defined in two— 
and three-dimensions by equations (3.39), (3.41) and (3.42) and the 
complementary functions are inferred from equation (3.7). The 
boundary stress calculation discussed in section 2.3.2 is used in 
place of the integral equation to represent the complementary function 
of stress for a point on the boundary. This avoids strong 
singularities in the integral equation associated with a boundary 
point . 

Selection of Parameters - The formulation presented thus far is 
relatively straight forward, however, the following areas demand 
special attention to guarantee high accuracy in the results. 
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The boundary mesh chosen for a specific problem should be fine 
enough to produce satisfactory results for elastostatic analysis. 
Fictitious function nodes should be introduced at all boundary 
collocation nodes and at points where interior quantities are to be 
evaluated since the interpolation of the global shape function is most 
accurate at these locations. Additional (fictitious function) nodes 
may be used through the interior of the body for a better 
representation of the particular integrals. A pattern of interior 
nodes consistent with the fineness of the boundary mesh is 
recommended, however, the number of additional points is dictated by 
the complexity of the temperature distribution. In any case, a 
uniform or neatly graded pattern is recommended for best results. 

It is obvious from equation (3.40) that the largest term in each 
column of matrix K falls on the main diagonal (when symmetric node 
ordering is used). The magnitude of the off diagonal terms will not, 
in general, be strongly banded unless special attention is paid in 
ordering the nodes according to their proximity. In any case, errors 
may show up when a large number of fictitious nodes are used. To 
eliminate this problem, matrix K can be optimized by the choice of the 
constants b n * These constants are dynamically calculated 
(independently for each n) to scale each column of matrix K so that 
the lowest value is forced to zero. This optimization and the use of 
double pivoting will reduce error in the solution of equation (3. S3), 
however, a price is paid for this matrix conditioning, since K will be 
rendered nonsymmetric. 

Still greater precision can be achieved by adding a constant term 
to the function in equation (3.40). This is accomplished by setting 
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one of the constants b n to zero in one specific column (here, for 
simplicity, we choose the last column; n = N). The benefits gained 
are seen when equation (3.40) is expanded out. For a field point x, 
we have 

T(x) = 2dC 


+ 3(l+d) 


b lPl0(?i) + 


+ b n P n 0U n ) + ... + b N _ x p n—1 ®(?n-1 ) 


+ 2d0(? N ) 


(3.55) 


where 

C - 0(?j) + 0(? 2 ) + ... + 0(? n ) + ... + 0(? N _ 1 ) 

P n = the euclidean distance between x and £ n » and 

2d ®^Ij) is the (self-adjusting) constant term 

For a given temperature distribution, C is constant (independent of 
x) , whereas the term in brackets depends on the spatial varying 
function p R , and therefore, the term in the brackets dictates the 
shape of the function. Note, the last term containing 0(? n ) (for n=N) 
is independent of coordinates and therefore is seen as a constant term 
that can adjust according to the value of ®(? N ). 

The values of the nodal densities 0(£ n ) must adjust to satisfy 
the temperature T(x m ) in N nodal equations. In the presence of the 
last term, C is seen as arbitrary, since the last term 0(?jj) can 
adjust to accommodate any value of C produced by the other N-l 0(£ n ) 
values. This allows the 0's for the terms in brackets to adjust to 
satisfy the temperature distribution in a free, relaxed manner. In 
the absence of the last term, C is no longer arbitrary and the 0(? n )'s 
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must take on values that force the term in brackets to balance C so 
that the temperature distribution is satisfied at nodal points. In 
essence, without the last terms, the global shape function may be 
forced to behave in a unnatural manner between the nodal points. 

3.5 EXAMPLES 

In this section three problems are presented to demonstrate the 
particular integral based, body force analysis. 

3.5.1 Gravity Stresses Around a Vertical Shaft 

The axisymmetric, gravity body force analysis is used to 

♦ 

determine the stress distribution in the soil near a vertical shaft in 
a halfspace due to its own weight. Eight quadratic elements (two for 
each side of a square region) with a total of sixteen nodes are used 
to model this geometry. In addition, stresses are calculated at four 
interior points. In figure 3.1, the BEM results are compared with the 
analytical solution (Terzaghi, 1943). Excellent agreement is 
obtained, even for this crude mesh. 

3.5.2 Stresses In a Rotating Sphere 

The axisymmetric centrifugal body force analysis is used to 
determine stresses in a solid sphere rotating about the vertical axis. 
Due to symmetry, only half of the body is modelled as shown in figure 
3.2. Four elements with a total of nine nodes are used. Note that 
the centerline is not discretized. The nodes along the r axis are 
constrained against displacement in the z direction. The only 
loading applied is that due to a rotation about the z axis. The 
BEM results are compared to theoretical and finite element 


90 



(Zienkiewicz, 1977) results in figure 3.3 and 3.4. Again, excellent 
agreement is obtained. 

3.5.3 Thick Cylinder Subjected to Thermal Load 

The axisymmetric, two- and three-dimensional thermoelastic 
analysis by particular integrals is tested on a thick cylinder under 
plane strain conditions. The boundary element discretizations for the 
problem are shown in figure 3.5. The thick cylinder, with an inner 
radius of 10 inches and an outer radius of 20 inches, is assumed to 
have the following material properties: 

E = 2600 lb/in 2 

v = 0.3 

a = 0.001 in/in/deg F 

The temperature distribution 

r 1 

T(r) = 1,000 — - 70 (deg F) 

is applied to the cylinder where r^ is the inner radius and r is the 
radial distance. The resulting BEM solutions for displacement and 
stress are given in Figures 3.6 and 3.7, respectively. These results 
are compared with the analytical solution by Boley and Weiner (1960), 
and once again, excellent agreement is observed. 

3.6 CONCLUDING REMARKS 

A new boundary element procedure, based on particular integrals, 
was introduced for the treatment of body forces. A range of problems 
was solved and shown to correlate favorably with existing results. 


91 



The particular integrals satisfy the inhomogeneous differential 
equation exactly, and therefore, the need for a volume integral or 
additional surface integrals is eliminated. The (homogeneous) 
boundary integral equation together with the particular integral 
represents an exact statement of the problem, and any error in the BEM 
solution is the result of approximations and errors introduced in the 
numerical implementation and solution. 

The general particular integral formulation introduced in this 
chapter can be applied in other BEM analyses which involve 
inhomogeneous differential equations. One such area, inelastic stress 
analysis, will be taken up in the next chapter. 
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Figure 3.1 

Stress Distribution near a Vertical Shaft due to Self-weight 
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R 


Figure 3.2 

Axisymmetric Mesh of a Solid Sphere 
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Figure 3.4 

Tangential Stress through a Rotating Sphere at Z=0 
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Figure 3.5a 

Two-dimensional Mesh of a Thick Cylinder (Diameter Ratio 


1 : 2 ) 



Figure 3.5b 

Three-dimensional Mesh of a Thick Cylinder (Diameter Ratio 1:2) 



Figure 3.5c 

Axisymmetry Mesh of a Thick Cylinder (Diameter Ratio 1:2) 



RRDIRL DISPLACEMENT (in) 



STRESS (ps i ) 
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Figure 3.7 

Stress through a Thermally Loaded Thick Cylinder (Plane Strain) 
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CHAPTER FOUR 


INELASTIC BOUNDARY ELEMENT FORMULATION 
BASED ON PARTICULAR INTEGRALS 


4.1 INTRODUCTION 

In the previous chapter, a general boundary element formulation 
based on particular integrals was presented and applied to 
gravitational, centrifugal and thermoelastic analysis. The next 
logical step is the extension of the technique to other areas of BEM 
analysis. In the present chapter, the formulation is extended to 
problems where the body forces are generated by the initial stress 
gradient; specifically elastoplastic analysis. 

Two- and three-dimensional particular integrals, as well as 
axisymmetric particular integrals are developed and implemented in a 
general purpose BEM computer code. The final BEM equation system 
produced by this formulation is similar in form to that generated by 
volume integration. This allows the iterative and variable stiffness 
algorithms developed in earlier works (Raveendra, 1984) to be 
incorporated with this new particular integral formulation without 
modification. Furthermore, the formulations developed in the next 
chapter for elastic and inelastic inhomogeneous media will employ the 
present procedure. 

4.2 PARTICULAR INTEGRALS FOR INITIAL STRESS BODY FORCES 

The governing differential equation for a body subjected to an 
initial stress body force is expressed in terms of displacement u A as 


101 



(4.1) 


(X+|l)U j,ji + * u i ,j j " ff ij,j 

where X and |i are Lam& constants. 

The Initial stress may be induced by a 'lack of fit' stress, 
temperature change, inhomogeneities, inelastic effects, or a 
combination of these. In this chapter a boundary element formulation 
using particular integrals is derived for the solution of the above 
equation, and from this formulation, the incremental formulation for 
plasticity is inferred. * * 

The total solution of the above equation can be decomposed into a 
particular integral and a complementary function: 


c p 

= u i + u£ 


- t c + t p 

Z L + H 

(4.2) 


a ij " °ij + °ij 

The boundary integral derived in Chapter 2 for elasticity 
represents the complementary function since it satisfies the 
homogeneous centerpart of equation (4.1). Hence, equations for the 
boundary system, displacement, and stress can be expressed in matrix 


form as 

Fu c = Gt c (4.3a) 

0 ° = G*t c - F*u c (4.3b) 

Using equation (4.2), these equations can be rewritten as 

Fu - Gt + [Fu p - Gt p ] (4.4a) 

0 = G®t - F®u + [F^uP - G*t p + ©Pi (4.4b) 
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The particular integral vectors u p , t p , and « p of the above equations 
are related to the initial stress by virtue of equation (4.1). Once 
explicit relations are known, these functions are evaluated at nodal 
points allowing equation (4.4) to be solved in the conventional manner 
for a set of well-posed boundary conditions. 

Effort is now turned toward deriving explicit expressions for the 
particular integrals. 


4.2.1 Two- and Three-Dimensional Particular Integrals 

The two- and three-dimensional particular integrals for 
displacement are related to the Galerkin vector via 



(1-v) p 

p i,kk 


— F 
2 p k,ki 


(4.5) 


where v is Poisson's ratio. 

Substituting this equation into equation (4.1) renders a 
relationship between the Galerkin vector and the initial stress 
function: 


F i,kkj j " (1-v) a ij,j (4 * 6) 

In subsequent steps of this derivation, it will be advantageous 
for the implicit expression of equations (4.5) and (4.6) to be related 
by a second order tensor, rather than a vector. Therefore, a tensor 
function h A j is introduced where 


ij,mmnn ij 


(4.7) 


Substitution of this equation into equation (4.6) and simplifying 
yields an expression for the Galerkin vector in terms of this new 
function: 
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F i “ (1-v) h ij,j 


Finally, substituting this expression into equation (4.5) yields the 
desired particular integral for displacement 


P 1 . ) y. 

u i = H il,lkk " 2p(l-v) lm.ilm 


(4.9) 


The particular integral for strain, stress and tractions are 


found using the following relations 


- I <u i.j + u j,i> 


_e p ° 
D ijkl 8 kl “ ®ij 


fc i ~ °ij n j 


(4.10a) 


(4.10b) 


(4.10c) 


where D e is the elastic constitutive relation given by equation 

J. J K i. 

(2.7a). It is important to note that the initial stress must be 
subtracted out of equation (4.10b) in order to produce the correct 
particular integral for stress. The complementary function for stress 

and strain, on the other hand, are related directly by 1 * e *' 

c _ _e c 
ij “ D ijkl 8 kl* 

Employing equations (4.10), the following expressions are 
obtained for the particular integral for strain and stress: 


8 ij = 2ji (h il,lkkj + h jl,lkki “ (1-v) ^m.ijlm^ 


(4.11) 


°ij = (1-v) h ml,mlkk 8 ij + h il,jlkk + h jl,ilkk 


(1-v) h lk,ijkl " h ij ,llkk 


( 4 . 12 ) 
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In passing we note the above relations for particular integral 
are derived in terms of initial stress, which is consistent with 
Chapter 2 . An initial strain e° ^ formulation is possible assuming 


h 


ij ,mmnn 


e ij 


(4.13) 


where the associate particular integral for displacement is 


u ? " TI^vT h kk,ijj + 2 h ij , jkk " TT^T h kj,ikj (4.14) 

The particular integrals given by equations (4.7), (4.9), (4.11) 
and (4.12) have little practical use in this implicit form. However, 
by applying the global shape function concept of the previous chapter, 
an explicit formulation can be developed for analysis of problems with 
general initial stress distributions. 


Global Shape Function - Tensor h^j(x) can be expressed in terms 
of a fictitious tensor density 0^j(£) as an infinite series using a 
suitable global shape function C(x,£): 


h nl <x) = I c(x -V (4 - 15) 

n=l 

m,l = 1,2 for two-dimensions 
m,l = 1,2,3 for three dimensions 

Several functions were considered, however, the best results were 
obtained with the following expression. 

C(x,? n ) . a£ [p 4 - b n p 5 ] (4.16) 

where 

A q is a characteristic length, 

p is the euclidean distance between the field point x and the 
source point ? n . 
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t> n is a constant, chosen in a manner described in section 
3.4.3. For the present discussion b R can assume the value 
of unity. 

All distances are non-dimensionalized by a characteristic length 

V 

The unknown fictitious densities are related to the initial 
stress through equation (4.7). Substituting equation (4.15) into 
equation (4.7) leads to 

CO 

°?m * ®lm<V <4 ' 17) 

n=l 

where 

K(x ' 5 n> = C ,mmnn (x ’S n > = a " bp 

a = 8d(d+2) , b = b n 15(d+3) (d+1) 

d = 3; for three-dimensional analysis, and 

d = 2; for two-dimensional (plane strain) analysis. 

The particular integral for displacement can be found by 
substituting equation (4.15) into equation (4.9). 


00 

“i<*> - 

n=l 

where 

^ml^'Sn 1 = A o [(Ci+<liP) <yiSi m + y m Sii) + <C2 +d 2P ) yi 8 im 

+ r y iyiy m ] 
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y l ■ 1*1 - <5„>il 


_ -8 
C 1 2|i(l-v) 


b 15 
n 

d l 2(i (1-v) 


c 


2 


C 1 


8(d+2) 
+ 



b 15(d+3) 
n 


The particular integral for strain and stress can be found by 
substituting equation (4.18) into equation (4.10). Alternatively, 
these functions can be found directly by a substitution of equation 
(4.15) into equations (4.11) and (4.12). Employing the latter method, 
the particular integral for stress can be written as 


°13 <"- 19 

where 

S ijlm (x 'V = (e 2 +f 2 p)8 jm 8 il + <e 3 +f 3 p) 8 ij 8 lm + (e 4 +f 4 p)8 im 8 jl 

+ V (y j y m 8 il + y i y j 8 lm +y i y m 8 jl ) 

^2 ^3 ^5 

+ T (y i y l 8 jm +y j y l 8 im ) + T y l y m 8 ij + J y i y j y l y m 


e^ = 2(iCj 
e 2 = ^( Cj + Cj ) 

= e 2 +X[c^(d+l)+c 2 l 
e 4 - e 2 ~ a 


f l “ “ f 5 “ 2pd l 
f 2 = ^(dj+dj) 

f 3 = fj + X[d 1 (d+2)+d 2 l 


Finally, the particular integrals for traction are found using 
equations (4.10c) and (4.19). 


107 



Once again the plane stress formulation is derived from the two- 
dimensional plane strain formulation by substituting the modified 
material constants defined in equations (3.16) and (3.17) into the 
plane strain equations. 

4.2.2 Axisysnetric Particular Integrals 

Although the methodology for axisymmetric BEM analysis by 
particular integrals is similar to that of two- and three-dimensions, 
the functions, of course, are different. In this section, an 
axisymmetric form of particular integrals is presented for 
axisymmetric initial stress body forces. 

t 

The axisymmetric particular integrals can be derived by one of 
the two methods described in Section 3.4.2. 

These methods are: 

1. A direct derivation utilizing the axisymmetric form of the 
differential operators in equations (4.7) and (4.9), and a 
suitable axisymmetric shape function. 

2. A derivation from the three-dimensional particular integral 
utilizing appropriate cylindrical coordinate and tensor 
transformations, and an analytic integration in the angular 
(0) direction. 

The second method has been adopted in this dissertation, and the 
resulting axisymmetric particular integrals for initial stress, 
displacement, and strain are given below. 
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< 


a ° r (x) 


(X) 

zz 


"rz<*> 


°©e< x) 


■i 




K U 0 0 K ?4 


k ; 2 o 


—r / 


0 


K « 0 0 ^ 


»rr'5„> 




0„<5 n > 

rz 




(4.20) 


'»?(*)’ 

00 

= \ 

U 11 

»; 2 

D ?3 

_n 

D 14 

< 

n=l 





uP(x) 

V * 


D n 

U 21 

D 22 

D 23 

D 24 


®rr^n* 


®zz^n* 


®rz<?n> 


e> ee (? n ) 


(4.21) 


£ r (x) 


P n 

E 11 

E 12 

E n 

e 13 

E n 

e 14 


»rr<V 

zz (x) 

00 

= \ 

E 21 

P^ 

e 22 

P n 

e 23 

F n 

e 24 


®zz^n^ 

rz<x) 

> n=l 

P n 

E 31 

p** 

E 32 

P n 

E 33 

E ?4 


> 

«W?n> 

ee (x) 


P n 

E 41 

P^ 

E 42 

F n 

E 43 

e ; 4 


0 ©O (5 n ) 


(4.22) 


where k”j, d”j and are subsequently defined. 

The particular integral expression for stress and traction are 

derived from equations (4.20) and (4.22) using equations (4.10b) and 

(4.10c). In the present implementation. *^j(x) is evaluated using 

equation (4.22) and the resulting values are substituted in equations 

(4.10b) and (4.10c) for the numerical determination of for the <r?,fv) 

ij' A ' 

and t^j(x). 
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The following notation is used in defining the axisymmetric 


particular integrals. 

r x = x r radial coordinate of field point 

z = v axial coordinate of field point 
X Z 

= U r ) n radlal coordinate of the fictitious density node i 

z § = ^z^n axda l coordinate of the fictitious density node n 

Note, the subscript n will be dropped for simplicity, 
coordinates are nondimensionalized by a characteristic length A q . 

z = z x - z 5 

R = [<r + r.) 2 + z 2 ] 172 

X S 

m = 4 w r2 

m 1 = m-1 

K = K(m) is the complete Elliptic Integral of the first kind 
E = E(m) is the complete Elliptic Integral of the second kind 
E and K are defined in Appendix III.D 

h l = 0 

h 2 = ji 

-1 

h 3 2|i( 1-v) 
a l = h 3 715 


All 
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a 2 “ h 3 /15 


a. = 


h 2 /3 + aj 


b. = 


h 3 /8 


b„ = 


h 3 /8 


b, = 


3h 2 /4 + bj 


b, = -b, 


General Form (r £ r. and/or z f z P ) 
£ . £ £ £. 

c i = K 

c 2 = (E-n^IO/m 

° 3 = [ 2E ( 2m-l ) + 1^(2-3111)] /3m 2 
= [4(2m-l)c 3 + SmjCjl/ 5111 
Cj = [6(2m-l)e^ + 5m^c 3 ]/7m 
Cg = I8(2m-l)c 3 + TmjC^l/Sm 

= [2(l+m 1 )E - mjKl/3 
e 2 = [4(1 + m 1 )e 1 - SmjEl/S 

~ E/nij 

d 2 = (K-E)/m 

d 3 = [Ed+fflj) - 21^] /m 2 
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2 

d 4 = [£(2-710! - 3mJ) + 

dg = [e 2 + 6m^E - 4m 1 e 1 - 4m 3 K + m* d^/m 4 

F 2 = 4Cj/R 

F 2 = -4(2c 2 - c x )/R 

F 3 = 4(4c 3 - 4c 2 + CjJ/H 

F 4 = -4(8c 4 - 12c 3 + 6c 2 - c-^/R 

G! = 4RE 

G 2 = 4mR (c 3 - c 2 ) 

G 3 = 2R [E-m (2c 4 - 3c 3 + c 2 )J 

Tj = 4d 2 /R 3 

3 

T 2 = -4(2d 2 - d^)/R 
T 3 = 4(4d 3 - 4d 2 + djJ/R 3 

3 

T 4 = -4(8d 4 - 12d 3 + 6d 2 - djJ/R 

T 5 - 4(16d 5 - 32d 4 + 24d $ - 8d 2 + d^/R 3 

If r x = 0 and/or r^ = 0 and z x t z^, then the following substitutions 

should be made for F . g and T 

a' a a 

F x = 2n/R 
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P 3 = V 2 



G x = 2nR 


G 


2 


= 0 


G 


3 


= JtR 


T 

T 

T 

T 

T 


= 2n/R 


= 0 


Tj/2 


= 0 


3T 3 /4 


Particular Integral for Initial Stress 
K n (x,5> = 7t - 3G 3 

K 14 (x,§) = jt - 3(0^) 

K 22 (x,lj) = 2n - 3G 2 

k 33 (x,^) = -3G 2 

K 41 (x,§) = Jt - SCG^Gj) 

* 44 (x,V = n — 3G 3 
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Particular Integral for Displacement: 

D 11 (x,§) = A q { n( 2a 2 +a 2 +a 3 )r x - Vr^-r^) 

- (b 2+ b 3 ) <r x G 3 -r ? G 2 > - b 2 tr x (r 2+ 2r 2 )F 3 - rj^F, 

- r ? (2r 2 *r 2 )F 2 - r x r 2 Fj] ) 

D 12 (x. 5) = A 0 { 2ita 2 r x - b 1 (r x 0 1 -r 5 G 2 ) - bjIr/j-i-jF,)* 2 ) 
D 13 (x,5) = -A 0 z C 2b 2 C(r ^ +r ^) F 2 - r x r 5 (F l +F 3 )] + (b 2 +b 3 )G 2 1 J 
d 14 (x,?) = A Q { 7 r( 2 a 1 +a 2 +a 3 )r x - (b 2 +b 3 ) (Gj-G 3 )r x 

- b i<r x Gi-r 5 G 2 ) - V* Ir^F^) - r ? (F 2 -F 4 )] J 

D 21 (x ' 5 ) = V { 2,ta l ■ b l G l " b 2 (r x F 3" 2 r x r 5 F 2 +r 5 F l ) } 
D 22 (x,?) = A^z { 2n(a 1 +a 2 +a 3 ) - (b 1 +b 2 +b 3 )G 1 - b^Fj^ } 

D 23 (x,^) = -A q { 2 n(a 2 +a 3 )r^ + (b 2 +b 3 ) ( r x G 2 -r £ G l^ 


+ 2b 2 2(r x f ' 2 - r ' ? F l ) 1 


D 24 (x,Jj) = A q z { Ina.^ - bjGj - b 2 r x (Fj-Fj) J 


Particular Integral for Strain: 

Eii<x.«> = n(2a 2 +a 2 +a 3 ) - bjGj - (b 2 +b 3 )G 3 

+ (3 V b 3 )r x r 5 F 4 - ' b l r 5 + b 2 <4r 2 +3r 2 > 4 b 3 (r 2 w 2 ) IF, 

+ V{ (2bj+5b 2 -b 3 )F 2 - Ibjrl+bjrllPj. 

+ b 2 Cr 2 r 2 (T 1 *T 5 ) - 2r x r ? <r 2 «- 2 HT 2 *T 4 > * <r^r 2 r 2 4r<)T 3 l 
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E 12 (x, 5> = 2naj - bjGj - b^F^ry^Fj+r^) - bjZ^ 

* V 2(r x T r 2 ‘y ? V'1V 

E 13 (x, 5> = (3b 2 + b 3 )r^zF 3 - ( 5b 2 +b 3 >r x zF 2 + 2b 2 zr^F 1 

* 2 V * ^l* 2r l )T 2 - '■x r ? T l 1 

= jt ( 22t^+a.2 + 2L3 ) ““ bjGj (^2*^3) (G^— Gj) 

' b l (r 'x F r 2r x r ? F 2 +r 5 F 3 ) - <b 2 +b 3> ,r x F r r x r 5 F 2- r x F 3 +r x r 5 F 4 ) 


- b 2 r x [ Sr^Fj-Fj) - 2r { (F 2 -F 4 > 1 


+ b 2 r x [ r x ( VV - 2r X*V W + ^'W 1 


E 21 (x,?) = 2 jta 2 - bjGj + b 2 z 2 (r^T 3 -2r x r 5 T 2 +r|T 1 ) 


- b l^ F l - b 2 <''x F 3- 2r x r ? F 2 +r 5 F l> 

E 22 (x,?) = 2n(a 1 +a 2 +a 3 ) - (b 1 +b 2 +b 3 )G 1 - [bj+4b 2 +b 3 ] z^Fj+bjZ^Tj 

E 23 (x.?) = -z(b 3+ 5b 2 ) <r x F 2 -r ? F 1 ) * 2b 2 z 3 (r x T 2 -r ? Tj) 

E 24 <x.{) - 2xa 3 - bjGj - bjZ^j - b 2 r 2 (F 2 -F 3 ) + b 2 z 2 r 2 <Tj-T 3 > 


E 3l(x,?) = 2 { “Z t (3b 2 +b 3 )r x F 3 “ (2b 1 +3b 2 +b 3 )r § F 2 +2b 1 r x F ]L ] 

* 2b 2 z t - r x r ?T 4 * r x (r 2 + 2r 2 )T 3 - r ? (2r 2 + r 2 )T 2 ♦ ry 2 ^ ] ) 


E 32 (x,$) = 2 C -z t(2b 1+ 3b 2 +b 3 ) (r x F 1 -r 5 F 2 )] + 2b 2 z 3 (r x T 1 -r ? T 2 ) } 
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E 33 (x '5> = \ { " 2(b 2 +b 3 )G 2 + (3b 2 +b 3 )[(r x r 5 (F 1 +F 3 ) - (r 2 +r^f-z 2 )F 2 l 
- 4b 2 z 2 [lyyTj+Tj) - <>'x +r | )T 2 1 } 

E 34 (x,§) = | C -2z ^(r/j-r^Fj) - z(3b 2 +b 3 )r x (F 1 -F 3 ) 

+ 2b 2 z r x tr x (T r T 3 ) - r ? (T 2 -T 4 )J } 

E 4^(x,£) = £) /^o r x 

E 42 (x,£) = D 12 (x,?)/A Q r x 
E 43 (x,?) = D 13 (x,£)/A Q r x 
E 44 (x,?) = D 14 (x,5)/A 0 r x 
If r = 0, then 

E 41<x»$> = E 11 (x,§) 

E 4 2<x,§) = E^ 2 (x»^) 

E 43 (x,£) = E 13 (x,^) 

E 44^ x '^ = E 14^ X ’^ 

Singular Form (r % = r £ and z x = z g ): 

The general form of the axisymmetrie particular integrals are 
singular when r x = and z x = z^. In this circumstance, the limiting 
form of these functions, given below, should be used. 
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G 3 = 1G^/1S 

7 1 - 4 
T 2 = -4/3 
T 3 = 28/15 

Particular Integral for Initial Stress: 

K 11 (x, 5) = n - 3G 3 

k 14 (x» 5> - * - StGj-Gj) 

K 22 (x,?) = 2 ji - 3Gj 

K 33 (x,5) = -3G 2 

K 41 (x,?) - Jt - SCG^) 

K 44 (x, 5> - it - 3G 3 

Particular Integral for Displacement: 

= A Q r x { 2na 1 + n(a 2 +a 3 > - b 1 (G 1 -G 2 ) 

- (b 2+ b 3 ) (G 3 -G 2 ) - b 2 r x (T r 2T 2+ T 3 ) ) 

D 12 (x,?) = A 0 r x { 2na 2 - b 1 (Gj-Gj) } 

D 13 (x,5) - 0 
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D 14 (x, 5) = A c r x { 2na 2 - bjtGj-Gj) ~ b 2 r x^ T l T 3* 

+ [(a 2 + a 3 )n - (bj+bj) (Gj-Gj)! ) 

D 21 (x,?) = 0 

d 22 ( x,5 > = 0 

D 23 (x,§) = -A 0 r x { 2«(a 2 +a 3 ) + (b 2 + b 3 ) (Gj-G^ } 

D 24 (x,?) =0 

Particular Integral for Strain: 

E n (x, 5 > = n(2aj+a 2 + a 3 ) - bjGj - (b 2 +b 3 )G 3 

+ r x t(3b 2 +b 3 ) (T 2 -T 3 ) - (b 1+ b 2 ) (W + \ b 2 (T^Tj+Tj)] 

E 12 (x, 5) = 2na 2 - b^ - b^tTj-Tj) 

E 13 (x,?) = 0 

E 14 (x,£), = Jt(2a 1 +a2 +a 3 > ~ b l G l “ * b 2 +b 3* * G 1 -G 3* 

" r x {b l (T l _T 2 ) + <1 b 2 +b 3 )(T l~ T 3 ) + b 2 [3(T l +T 2 ) " 2<V T 3 )1} 

E 21 (x, 5) = 2na 2 - bjGj - b 2 r x (T 1 -T 2 ) 

E 22 (x, 5) = 2n(aj+a2+a 3 ) - (b 1 +b 2 +b 3 )G 1 

e 23 (x,|) = 0 

®24^ x ’^ = ^ na l ~ b l G l ~ b 2 r x^l + ^2^ 

E 31 (x,?) = 0 
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E 32^ x '^ 

= 0 

E 33 (x,£) 

- - < b 2 +b 3 ,0 2 * <f b 2 + 5 WW 

E 34 (x,5) 

- 0 

e 41 <x,5) 

= Djj(x>§)/A 0 r x 

^ 4 2^ x » 5) 

= D 12 (x,5)/A Q r x 

E 43 ( x,5) 

= 0 

E 44<x,5) 

= D 14 (x,£)/A 0 r x 

If r =0, then 

X 


e 41 <x,S) 

= E n (x,5) 

e 4 2 ( x.5) 

= E 12 (x * 5 ) 

E 43 ( x,5) 

= 0 

e 44 (x,5) 

= E 14^ x '^ 


4.3 NUMERICAL IMPLEMENTATION 

The present formulation is implemented in a manner analogous to 
Chapter 3. Essentially, the procedure consists of evaluating the 
relevant particular integrals at nodal points and solving equation 
(4.4) for these values and a set of appropriate boundary conditions. 
The particular integrals are a function of initial stress. However, 
the initial stress of an inelastic analysis is unknown and must be 
determined as part of the solution process. Therefore, it is 
necessary to assemble the equation system in a manner that will admit 
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to an inelastic solution algorithm. Both the iterative and the 
variable stiffness algorithms of Chapter 2 can be employed without 
modification if the assembly process presented below is utilized. 

For the purpose of numerical evaluation, the particular integral 
series solutions of the previous section are truncated to a finite 
number of terms. Note, the particular integrals of each region are 
evaluated independent of the other regions, and calculated only in 
regions where inelastic effects are anticipated. 

Equations for u p and t| are written for the boundary nodets and 
are expressed in matrix form for each region as 

u p = D 0 

(4.23) 

t p = T 0 

in which D and T are matrices of the order (f*M) by (g*N) where M is 
* 

the number of boundary nodes, N is the number of terms in the series, 
f is the number o£. degrees of freedom of the analysis, and g is the 
number of independent stress components . 

The initial stress, evaluated at the N nodal points, is expressed 

as 

a° = K 0 <4.24) 


in which K is a well conditioned (g*N) by (g*N) matrix. Post- 
multiplying equations (4.24) by K -1 yields 

0 = K" 1 o° (4-25) 

Back substituting this equation into equation (4.23) renders 


u p - D K" 1 a° 
t p = T r 1 a° 


(4.2 6) 
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Similar particular integral expressions can be written for 
displacement and stress at interior points of interest. The 
particular integral solution for stress, corresponding to nodes of 
equation (4.4b), is written in matrix form as 

o 1 * = S K" 1 a° (4.27) 

Substituting equations (4.26) and (4.27) into equation (4.4) and 
rearranging yields 

Gt - Fu + Bo° = 0 

(4.28) 

a = G®t - F°u + B V* 

where 

B = -[GT - FD]K -1 
B® = -[G°t - F°D] K -1 


The equations of (4.28) are in the same form as those obtained by 
conventional volume integration. After the above equations are 
generated for all regions in a problem, they are assembled in a manner 
described in Chapter 2. The final system equations are expressed as 


A b x = B b y + C b a° 


<r = k a x + B a y + C®«° 


(4.29) 


A solution for the inelastic, incremental form of this equation is 
obtained using one of the two algorithms presented in Chapter 2. 

Finally we note, two time-saving features that are employed in 
the present implementation. 
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First, careful observation of equation (4.15) reveals that the K 

2 

matrix (equation 4.25), in two- and three-dimensions, consists of N 
diagonal matrix blocks. All diagonal terms of a particular diagonal 
matrix block are the same value. By reducing each diagonal matrix 
block to a single term, the overall reduced matrix is inverted at a 
fraction of the cost of the whole. Expanding each term of the 
inverted matrix back to a diagonal matrix block produces the correct 
K -1 . 

Second, the stress for a boundary point is directly (and 
efficiently) calculated by the method described in section 2.3.2, 
instead of using the formal particular integral procedure of this 
chapter. 

4.4 EXAMPLES 

In this section, three examples are presented to demonstrate the 
particular integral based, inelastic analysis. 

4.4.1 Three-Dimensional Analysis of a Cube with Hardening 

The three-dimensional, inelastic particular integrals are tested 
on a unit cube in tension with plastic strain hardening. The material 
constants are (in consistent units): E = 100, v = 0.3 and h = 50.0. 

The surface of the cube is discretized (figure 4.1) using six (eight- 
noded) boundary elements. The particular integrals are defined using 
the twenty boundary nodes and one interior node located at the center 
of the cube. The displacement in the axial direction on a face 
opposite the fixed end is shown in figure 4.2 for increasing tension. 
The particular integral based results are in excellent agreement with 
the analytical solution. 
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4.4.2 Axisymmetric Analysis of a Thick Cylinder 

A thick cylinder (1:2 ratio) subjected to interior pressure is 
analyzed using the axisymmetric, inelastic particular integral 
formulation under the plane strain condition. The axisymmetric mesh, 
shown in figure 4.3, has ten quadratic boundary elements. Twenty- 
three nodes are used to define the particular integral domain 
representation. The load-displacement response at the outer surface 
is in good agreement with the analytical solution (Hill, 1950) as 
shown in figure 4.4. The hoop strain at the inner and outer surface 
is shown in figure 4.5 for increasing pressure. Once again, the 
results are in good agreement with the analytical solution. 

4.5 CONCLUDING REMARKS 

The particular integral based boundary element formulation 
introduced in Chapter 3 was successfully extended to inelastic 
analysis in the present chapter. The system matrices produced by this 
method were assembled in the same form as those created using the 
volume integral based method, and therefore, the two inelastic 
solution algorithms, presented in Chapter 2, could be applied without 
any modifications. 

The method was demonstrated for a number of problems and 
excellent agreement was obtained in comparison with existing results. 
Practical applications of the present analysis will be presented in 
Chapter 6. 
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Figure 4.1 

Three-dimensional Mesh of a Cube 
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Figure 4.2 

End Displacement of a Cube (with Strain Hardening) under Tension 




Figure 4.3 

Axisymmetric Mesh of a Thick Cylinder for a Particular Integral based 
Plasticity Analysis (Diameter Ratio 1:2) 
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Figure 4.4 

Radial Displacement of the Outer Surface of a Thick Cylinder under 
Internal Pressure (Plane Strain) 
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Figure 4.5 

Hoop Strain in a Thick Cylinder (Plane Stress) 
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CHAPTER FITE 


BOUNDARY ELEMENT FORMULATION FOR INHOMOGENEOUS MEDIA 
5.1 INTRODUCTION 

One deficiency of the boundary element method is its inability to 
readily handle continuous material inhomogeneities. However, a method 
analogous to the plasticity procedure can be employed to overcome this 
difficulty. In this procedure, the effects of the inhomogeneities in 
the material are incorporated in the boundary element system through 
an initial stress body force. 

In this chapter, formulations for elastic and inelastic analysis 
of inhomogeneous media are presented. The inhomogeneities may occur 
as a result of spatial variation in material parameters, or may be 
thermally induced when material parameters are assumed to be 
temperature dependent. 

Successful inhomogeneous BEM formulations have been implemented 
by Butterfield (197 8) for potential flow problems and by Ghosh and 
Mukherjee (1984) for thermoelastic bodies. Both of these linear 
analyses were based on iterative procedures. The present analysis is 
a direct implementation which incorporates the inhomogeneous effects 
into the system without iteration. This leads to a direct solution 
procedure for elastic inhomogeneous problems and a single iteration 
process for plasticity analysis. The formulations are implemented in 
a two-dimensional and axisymmetric, multiregion, general purpose 
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system. Since the procedure involves an initial stress body force, 
the volume must be appropriately represented in both elastic and 
inelastic analysis. Either the conventional volume integral or the 
new particular integral representation may be used for this purpose. 

5.2 BEM FORMULATION FOR ELASTIC INHOMOGENEOUS MEDIA 

In Chapter 2, the boundary integral equations are derived via the 
Betti reciprocal work theorem. Although, this theorem is valid for 
continuums of general material, the adoption of the Kelvin point force 
solution renders the resulting integral equations valid only in 
homogeneous, isotropic material. Nevertheless, material 
inhomogeneities can be incorporated in the standard boundary integral 
equations through a volume integral or the particular integral 
equivalent . 

The differential equation expressed in terms of stress (equation 
2.1) is independent of material, and therefore, valid for 
inhomogeneous material. 

*ij.j + f i - 0 <*•*> 

The strain in an elastic, inhomogeneous, isotropic body is related to 
stress via 

ij " D ijkl c kl 
or 

m _ e,L 
e kl “ c ijkl a ij 

where : 

the superscript L indicates the use of a (local) constitutive 
relation that varies with position, i.e.. 


(5.2a) 

(5.2b) 
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D ijkl - >•<*> 5 ij 6 kl * 8 ik S J1- and 

mechanical strain e®^ is defined in terms of total strain 


8 kl " 8 kl " 


8 kl o(x) T 


(5.3) 


Substituting the strain-displacement equation (2.4) and equation 
(5.2) into equation (5.1) yields the displacement form of the 
equilibrium equation. Since material parameters are a function of 
position, the resulting expression will include spatial derivitives of 
these parameters. These derivitives can be avoided, however, by 
recasting the constitutive relation in an alternate form as described 
below: 

The mechanical strain e®^ is divided into two parts; a 
homogeneous strain and an inhomogeneous strain j . 


e 


m h 
ij = e ij 



(5.4) 


The homogeneous strain is defined with respect to global reference 
material parameters X G and n G that do not vary with position. 


e h - r e,G n 
e ij " C ijkl "kl 

or 

_ _.e,G h 
ff ij “ D ijkl 8 kl 


(5.5) 


(5.6) 


where cVHi, are the global compliance tensors. The 

ljKi ijKX 

inhomogeneous strain is simply the supplement part 


p i _ m h 
e ij " 8 ij " 8 ij 


(5.7) 
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Using equation (5.4) in equation (5.6) yields an expression for stress 
in terms of the global compliance tensor: 


ij 


= D 


e,G m 
ijkl 8 kl 


i 

a kl 


(5.8) 


where 



= D 


e,G i 
ijkl e kl 


(5.9) 


or using the definition of mechanical strain, a relation between 
stress and total strain can be expressed as 


ff ij " D ijkl e 


kl 


t 

“ 0 kl 


i 

^kl 


where 


(5.10) 




8 kl T 


(5.11) 


Note in this calculation the local value o(x) is assumed, but the 
constitutive relation is with respect to global reference parameters . 

Substitution of equation (5.10) and the strain-displacement 
relation into equation (5.1) leads to the displacement equilibrium 
equation expressed in terms global material parameters: 




U j >ij 




(5.12) 


where 


«° t' i 
ij = ®ij + c ij 

c c 

Note, since X u and n are not a function of position, they are removed 
from within the differential operators. 
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Equation (5.12) is the extended form of the Navier equation, 
expressed in terms of the global reference parameters. The boundary 
integral equations satisfying equation (5.12) were presented in 
Chapter 2. The expression for displacement (neglecting body forces) 
is repeated here. 

C ij (§)u i (?) = j [G ij (x,§)t i (x) - F iJ (x,§)u i (x)l dS(x) 

S 

+ i B ikj^ x »$) a ik^ dV * x ) (5.13) 

V 

where the kernel functions assume the global reference parameters. 
Similar expressions for stress and strain can be written. 
Alternatively, a solution satisfying equation (5.12) also can be 
formulated using the particular integral procedure presented in the 
previous chapter. Regardless of the formulation employed, the final 
form of the system equations are the same for both the volume integral 
and the particular integral based procedures. The solution to the 

boundary value problem requires the knowledge of the inhomogeneous 
part of the initial stress. To this end, equation (5.2) is set equal 
to equation (5.8), and upon rearranging, the following relation is 
obtained . 



.e , L * 
5 ijkr 


m 

e kl 


(5.14) 


This expression is analogous to equation (2.52c) of the 
plasticity formulation and it plays a similar role in the iterative 
procedure for inhomogeneous problem. However, unlike plasticity, the 
elastic, inhomogeneous analysis is not incremental. 


134 



The inhomogeneous problem can be solved using the following 
iterative procedure. The integral equations for the boundary system 
and interior stresses are integrated and assembled in the usual manner 
(Chapters 2 and 4). is initially assumed to be zero. The 

boundary problem is solved under this condition and the resulting 
stresses (and strains) are determined at cell nodes (or particular 
integral nodes). Using equation (5.14)» the inhomogeneous part of the 
initial stresses are determined at these same nodes. In addition to 
the external loading, the calculated values for the inhomogeneous 
(initial) stress are applied to the system. The problem is resolved 
and new values for are determined. This process is repeated until 
a convergenced solution is obtained. 

The rate of convergence is dependent on the degree of 
inhomogeneity, the choice of global reference parameters, and the 
geometry and loading of the problem. In some cases the convergence 
may be slow or may not converge at all. For this reason, a direct 
algorithm, analogous to the variable stiffness plasticity method of 
Chapter 2, is developed. In this procedure, the unknown, 
inhomogeneous part of the initial stress is eliminated from the 
boundary system permitting a direct solution. 

The present formulation is implemented on a nodal basis, and 
therefore, the subsequent derivation is carried out on the discrete 
nodal equations rather than on the (continuous) integral equation. 
First, the integral equations are integrated for a system of nodal 
equations and assembled in a manner described in Chapter 2 or Chapter 
4. Since the assembled form of the boundary equations are similar for 
both the volume integral base method and the particular integral 
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formulation, no distinction is necessary in the present derivation. 
These assembled equations are repeated here in matrix form: 

A b x = B b y + CV + C b o*- (5.15a) 

<r = k°x + B°y + cV + CV (5.15b) 

where a stress equation is written for every node where the unknown 

exist. Therefore, C° is a square matrix. 

Substituting equation (5.2b) into equation (5.14) we arrive at 
the nodal relation 



(D 


e,G _ n e,L . _e,L 
ijkl D ijkl' °klmn a mn 


(5.16) 


In the present work, we will assume the modulus of elasticity may 
vary with position, but Poisson's ratio is constant over the domain. 
This is a realistic assumption for most materials. Upon invoking this 
assumption, the following simplification is possible. 


(D 


e , G 
ijkl 


- n e * L } 
D ijkl> 


-,e,L 

■'klmn 


E G -E L 


6 im 6 jn 


(5.17) 


where E® and are global and local moduli of elasticity, 
respectively* Substituting this result in equation (5*16) and 
rearranging yields 


c ij 



(5.18a) 


where 


H = 


e g -e l 


(5.18b) 
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Applying this relation to the nodes of equation (5.15), a matrix 
expression can be written: 

o = Ebo 1 (5.19) 

where I represents the identify matrix and h represents a nodal vector 
of corresponding nodal values of H. 

Substituting equation (5.19) into equation (5.15b) yields 

Iho* = A°x + B°y + C °a t + C®* 1 (5.20) 

and upon rearranging equations (5.15a) and (5.20) we have 

A b x = b b + cV- (5.21a) 

Do 1 = A°x + b ff (5.21b) 

where 

b b = B b y + 
b a = B w y + 

D = Ih - C a 

Equation (5.21b) can then be recast as: 

e 1 = D -1 (A°x + b a ) (5.22) 

Substituting the above equation into equation (5.21a) results in the final 
system equation: 

A*x - b* (5.23) 

where 

A* = A b - EA®, 
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b* = b^ + Eb ff , and 

B = cV " 1 

This equation can be solved for the unknown vector x using 
standard numerical techniques. Once x is found, equation (5.22) is 
used to calculate the inhomogeneous part of the initial stress, and 
equation (5.19) is used to evaluate the corresponding stress at these 
nodes. With vector x and o* determined, the displacement, stress or 
strain can be found at any point. 

In a multi-region system, the matrices of the above equations are 
block-banded. When the above formulation is efficiently implemented, 
this property leads to a reduction of work in the matrix 
multiplication. Furthermore, the inversion of matrix D is carried out 
one region at a time in a substructured analysis. 


5.3 BEM FORMULATION FOR INELASTIC INHOMOGENEOUS MEDIA 

The total inelastic strain increment is comprised of an elastic, 
plastic and thermal strain rate components: 


_ e . p . t 

c ij " 8 ij + e i-j + c ij 


(5.24) 


The elastic strain rate can be unceremoniously divided into two parts: 

*h *i 

a homogeneous strain rate e^; and a inhomogeneous strain rate e^. 


e ij = 8 ij + 8 ij 


(5.25) 


using global material constants, can arbitrarily be defined as 

8 ij = D ijkl ff kl (5.26a) 
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or 


(5.26b) 


'ij 


= D 


e,G h 
ijkl e kl 


i e h 

and 6^.^ is defined as the difference between and e!y 

In light of equation (5.25), equation (5.24) can be rewritten as 


»h _ tip 

8 ij " 8 ij ” 8 ij " 8 ij c ij 


(5.27) 


and substituting this result in equation (5.26b) leads to the 
following relation for stress rates. 


_ _ n e,G n e,G t _e,G i ^e.G p 

ij " D ijkl c kl “ D ijkl 8 kl " D ijkl 8 kl " D ijkl e ij 


(5.28) 


Rewriting the above equation in simpler form yields a global 
constitutive relation for inhomogeneous plastic media. 


'ij 


n e,G o 

D ijkl e kl _ ®ij 


where 


- 

V 

; ij * 


’ij 


+ a 


ij + ®ij 


) e »® 

; ijkl e kl 

e.G *i 
J ijkl c kl 


rP. = D e 'G 


ij 


ijkl e kl 


(5.29) 

(5.30a) 

(5.30b) 

(5.30c) 

(5.30d) 


The displacement rate formulation of the equilibrium equation for 
inhomogeneous media is derived by substituting the above constitutive 
relation and the strain-displacement relation into the stress rate 
equilibrium equation. The resulting equilibrium equation is 
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where X G ,ji G are the global material constants, and the initial stress 
rates are defined by equation (5.30). 

Once again, it is apparent that the boundary integral equation 
derived in Chapter 2 or the alternate particular integral formulation 
of Chapter 4, satisfies the above equilibrium equation. The thermal 
(initial) stress rates are calculated from the known temperature 
distribution, however, the initial stress rates for the plastic and 
inhomogeneous effects, are unknown and must be determined as part of 
the solution procedure. Therefore, it is necessary to derive a 
relation between these initial stress rates and the current rate of 
stress . 

A relation for the inhomogeneous (initial) stress will be 
derived first. In both the elastic and plastic regions, the elastic 
strain rate is related to the stress rate via the local elastic 
constitutive relation. 


or 


"id 


= D e '- L - 


ijkl e kl 


ij 


D ijkl (8 kl 


- e 


kl 


(5.32a) 

(5.32b) 


where 

• • 

8 kl = e kl ~ c kl ' and 

ej^ =0 in the elastic region 

(This equation is analogous to equation (5.2) of the previous 
section.) Equating this with equation (5.28) and rearranging leads to 



the desired relation 


or 


’ij 


U 


/ T\® > ^ —6 > L t ^ p \ 

(D ijkl " D ijkl ) (8 kl ” e kl ) 


ve,G 

J ijkl 


ve,L 

'ijkl' 


i U f^6 > b v / ID i 

(D4 41,1 — D4 4 t,T ) 


for the elastic region 


(5.33a) 

(5.33b) 


Next a relation involving the plastic (initial) stress is derived 
using the local elastoplastic constitutive relation (equation (2.50)): 


_ _ _.ep,L m 

ij ~ D ijkl e kl 


(5.34) 


Equating this to equation (5.29) yields 


• . • p - ♦ 

ff ij + °ij “ ^ D ijkl ” D ijkP e kl for the P lasfcic region. (5.35) 

The equations above suggest the use of a single iterative 
procedure in which the effects of the inhomogeneity and plasticity are 
treated together. In an iterative algorithm (similar to the one in 
Chapter 2) a small load increment is applied to the system, and the 
resulting elastic stress state is determined. Theoretically, this 
stress state is checked at nodal points through the domain and 
depending whether or not the yield condition is violated, either 
equation (5.33b) or equation (5.35) is employed. However, when 
inhomogeneities exist in the body, the boundary integral equations 
will yield a fictitious state of stress that is dependent on the 
global reference parameters, and therefore, the choice of equation 
(5.33b) or (5.3 5) is not clear. The proper choice, however, is 
critical since plastic deformation is irrecoverable and the use of the 
wrong equation will ultimately result in error. So, before the effect 
of plasticity can be accounted for, the actual inhomogeneous elastic 
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stress state must be determined by its own initial stress correction. 
Hence, a double iterative procedure is necessary to account first, for 
the inhomogeneity effect, and second, the plasticity effect. More 
specifically, between each iterative step of the plasticity algorithm, 
it is implied that a converge inhomogeneous state is achieved through 
a separate iteration. 

This double iteration is very costly and a scheme for 
accelerating the convergence (similar to that of Chapter 2) would be 
essential. Also, further study may reveal that additional savings can 
be obtained at later iterations, i.e., during initial iterations of 
each load step, a double iteration process is necessary, however, at 
later iterations, when the plastic regime is clearly defined, a single 
iterative procedure may prove viable. 

Nevertheless, a semi-direct method, based on the direct 
formulation of the previous section, is possible. In this procedure, 
the inhomogeneous effects are directly incorporated in the boundary 
system through a back— substitution of the stress rate equations, and 
the plastic deformation is incorporated in the system through an 
iteration process. Hence, the method entails only a single iteration 
process in which inhomogeneity is satisfied at all times. In the 
present work this approach is adopted. Therefore, it is necessary to 
separate the plasticity and inhomogeneities effects of equation 
(5.35). This could be accomplished by subtracting equation (5.33a) 
from equation (5.35), resulting in 


ur , = 


U 


(D ijkl " D ijkl )e kl 


+ (D 


e,G 

ijkl 


.e,L 

5 ijkl' 




D e,G 

ijkl 


*1 


(5.36) 
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This equation, however, is not an explicit relation between the 
plastic (initial) stress rate and the current rate of stress (or 
strain). Therefore, this equation is reformulated in an alternate 
manner consistent with Chapter 2. From equation (2.52c), the 
mechanical strain at a point in which yielding occurs is related to 
the plastic (initial) stress rate via 



/ n e »L # Lt id 

(U ijkl D ijkl ,8 ijkl 


(5.37) 


The bar indicates that the plastic stress rate is related to the 
plastic strain rate through the local constitutive relations as 
follows 


»P _ p6» L , n e,L _ — ep,Lv m 
e ij ^ijkl (D klmn D klmn' e mn 


(5.38) 


Substitution of this result into equation (5.30d) yields the desired 
form of the relation: 


gP in® » (* <-i®»L .*ep,L. m 

ij u ijkl L klmn ' D mnop D mnop ,e op 


(5.39) 


(This is a restatement of equation (5.36)). With equations (5.33) and 
(5.39), the semi-direct BEM procedure for inelastic, inhomogeneous 
media can be formulated. 

We begin with the system equations (from either Chapter 2 or 4) 


expressed in matrix form: 



Once again, the Poisson ratio is assumed not to vary spatially 
within a substructured region. Therefore, the following relation 
between the stress rate and the inhomogeneous (initial) stress rate 
can be written: 


*i 1 ‘ 
ff iJ = H a ij 


(5.41) 


where 


E g -E l 


or expressed in matrix form for the nodes of equation (4.40b) 


- Ih or 1 


(5.42) 


Substituting this into equation (4.40b) leads to 


= A°x + B a e b + C ff a b + 


(5.43) 


and rearranging equations (5.40a) and (5.43) yields. 


A b x = b b + C b ^- + C b > 


(5.44a) 


c 1 = D -1 (A a x + b a + C®oP) 


(5.44b) 


where 


b b «■ B b y + chffh 


b® = B a j + 


D = Ih - C c 
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Substituting equation (5.44b) into (5.44a) results in the final system 
equation: 

A*x = b* + C* (5.45) 

where 

A* = A b - EA ff 
b - b b + Eb* 

C* = C b + EC* 

E = 

This equation system for inhomogeneous inelastic material takes 
on the same form as the homogeneous, inelastic equation (see eq. 
2.54), and can be solved in a similar fashion using the iterative 
plasticity method. The stress at cell nodes (or particular integral 
nodes) can be found using equation (5.42) where the corresponding 
inhomogeneous (initial) stress o* is defined by equation (5.44b). 
After the unknowns on the boundary have been determined, the 
displacement and stress at any point in the body can be found with the 
appropriate discretized integral equation. 

A few comments are in order. First, when the inhomogeneities are 
caused by a prescribed variation in material parameters, the 
inhomogeneities remain constant and the boundary system need only be 
constructed once. However, thermally induced inhomogeneities vary 
with change of temperature and the boundary system should, in 
principle, be reconstructed at every load increment. Thermally 
induced inhomogeneities are, however, very mild and reconstruction is 
not necessary at every load step. Instead, reconstruction can be 
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delayed until the resulting error surpasses a prescribed tolerance. 

This brings us to the second point; efficiency. Although it is 
necessary to explicitly construct the A* matrix (and in the process 
construct the D -3 and E matrices), other matrix manipulations 
associated with the right hand side of equation (5.45) are unnecessary 
and can be handled more efficiently through a series of matrix/vector 
multiplications. More specifically, equation (5.45) should be set up 
in the following manner. 

A*x = b 1 + b 2 (5.46) 

where 

b 1 = b b + C b > 
b 2 = Eb 3 
b 3 = b a + CV* 

Note vector b 3 can be used a second time during the calculation of the 
inhomogeneous (initial) stress rates (equation (5.44b)) and subsequent 
calculation of the real stress rate (equation (5.42)). 

Finally, it is possible to formulate an entirely direct approach, 
for the analysis of inelastic, inhomogeneous media, by employing the 
variable stiffness plasticity technique together with the direct 
formulation for inhomogeneities. However, this approach is very 
costly and inefficient due to the large amount of matrix manipulations 
involved. For instance, the matrix multiplications in equation 
(5.44b) and (5.42) has to be carried out explicitly in order to form 
the inhomogeneous stress rate equation. And, in addition to the 
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construction of matrix k* of equation (5.45), the C* matrix also has 
to be explicitly constructed. Finally, once the inhomogeneous 
operations are complete, the matrix operations associated with the 
variable stiffness method must still be carried out. Therefore, the 
overhead involved in the setup of the completely direct procedure will 
inhibit its chance of becoming a viable alternative. This is 
particularly true for thermally induced inhomogeneities since it 
requires repeated construction of the inhomogeneous system. 

5.4 EXAMPLES 

In this section, three examples are presented to demonstrate the 
inhomogeneous formulations. Inhomogeneous variation is usually 
relatively mild in most problems. However, for purpose of 
demonstration, severe variations are assumed in these examples. 

5.4.1 Elastic Rod with Spatially Varying Modulus 

An elastic cylindrical rod shown in figure 5.1a has an elastic 
modulus that varies through the length. A uniform traction is applied 
to top, and the bottom surface is fixed in the axial direction. An 
axisymmetric discretization, shown in figure 5.1b, has three boundary 
elements and one volume cell. The material properties of the rod are 
(in consistent units): 


E(z) = 


100 

1+z 


= 0.3 


The global reference modulus is assumed to be unity E® = 1.0. The 
axial strain through the cube is shown in figure 5.2 for an applied 
traction of t = 100. Good agreement is obtained with the analytical 
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solution e = (l+z) . 
z 

5.4.2 Elastic Cube with Temperature Dependent Modulus 
An elastic cube with a temperature dependent elastic modulus is 
analyzed under plane strain conditions. For a cube in plane strain 
with boundary conditions shown in 5.3a and subjected only to thermal 
loading, the following relations are obtained 


e 


x 


(1+v) 

(1-v) 


aT 


c y 8 z e xy 0 


°x - ff xy ~ 0 


_ aT 

°y ~ a z ~ ~ E (1-v) 


In most materials, the modulus is usually assumed to vary as a 
linear function of temperature within a prescribed temperature range. 
In terms of temperature T this can be expressed as 


E(t) = E q + k(T-T 0 > 

where E q is the elastic modulus at the reference temperature T Q . The 
change in modulus in most real materials is usually quite small, and 
therefore, k is usually a very small value. However, in this example 
a severe variation of elastic modulus (k=l) is assumed for 
demonstration purposes. The assumed temperature distribution and 
material properties are (in consistent units): 


T(x) = lOOx 
a = 0.7 x 10" 3 
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V 


= 0.3 


E = 100.0 at T = 0.0 
o o 

k = 1.0 

Therefore, 

E(T) = 100 + T 
or 

E = 100 (1+x) 

Substituting these values in the expression for stress and strain 
yield 

e x = 0.13x 

<r y = -10x(l+x) 

and integrating the strain across the length yields an expression for 
displacement : 

u x = 0.065X 2 

The discretization of the cube, shown in figure 5.3b, has four 
boundary elements and nine particular integral nodes. The 
displacement and stress profile obtained from the analysis is shown in 
figures 5.4 and 5.5, respectively. Good agreement is obtained in 
comparison to the analytical solutions. 

5.4.3 Plastic Analysis of a Cube with Spatially Varying Modulus 
The plastic deformation of an inhomogeneous cube (figure 5.6a) in 
tension is analyzed under plane strain conditions. The 
discretization, shown in figure 5.6b, consists of four boundary 
elements and one volume cell. The material properties, described 
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below, assume a modulus of elasticity that varies as a function of x. 
Given in consistent units: 

E = 100/ (1+x) 

v =0.3 

$ 

o Q = 100.0 

h = 50.0 

* % 

E° = 1.0 

In figure 5.7, the displacement rate profile through 'the cube in 
the x direction (at y = 0) is given for loads: t x = 100. (yield), 

t = 110., t = 120., and t„ = 130. Good agreement is found in 

X. X X 

comparison with the analytical solution. A similar graph is given for 
the lateral displacement rate (at y = 0.5) in figure 5.8. Once again, 
good agreement, with the analytical solution, is obtained. 

5.5 CONCLUDING REMARKS 

An axisymmetric and two-dimensional, boundary element 
formulations were derived for elastic and inelastic, thermal, 
inhomogeneous media. The elastic formulation employed a direct 
solution procedure. The inelastic formulation used a semi-direct 
procedure based on the iterative plasticity algorithm presented in 
Chapter 2. The results obtained by this method for the example 
problems were in excellent agreement with analytical solutions. 

The present formulation can be extended to three-dimensional 
analysis without difficulty. 
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Figure 5.1a Figure 5.1b 

Inhomogeneous Cylindrical Rod in Axial Tension Axisymmetric Mesh of a Cylindrical Rod 

with One Volume Cell 
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Figure 5.2 

Axial Strain through an Inhomogeneous Cylindrical Rod in Tension 
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Figure 5.3a 

A Cube of Thermal-sensitive Material (Plane Strain) 


Figure 5.3b 

Discretiztion of a 2-D Cube using a 
Particular Integral Domain Reprsentation 
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LATERAL STRESS 
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Figure 5.5 

Thermally Induced Inhomogeneous Lateral Stress {a y and a z ) 
in a Cube (Plane Strain) 
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Figure 5.0b 

Discretization of a 2~D Cube with One Volume Cell 


Figure 5.6a 

Inhomogeneous Cube in Axial Tension 
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Figure 5.8 

Inhomogeneous Plastic (Lateral) Displacement in a Cube in Tension 
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CHAPTER SIX 


ADVANCED ENGINEERING APPLICATIONS 


6.1 INTRODUCTION 

In previous chapters, elementary examples such as cubes, spheres 
and cylinders are used to verify the accuracy and convergence of 
proposed formulations. These problems were chosen since their 
simplicity gave way to analytical solutions. 

In this chapter, more challenging problems of practical interest 
are considered. 

6.2 ELASTIC ANALYSIS 

6.2.1 Bending of a Circular Plate 

A clamped, circular plate subjected to a uniformly distributed, 
circular, patch load is shown in figure 6.1a. The plate is modeled by 
a five region, axisymmetric mesh shown in figure 6.1b. The modulus of 
elasticity is E = 210,000, the Poisson's ration is v = 0.3, and a/b = 
10 . 

In figures 2 and 3, the radial, tangential, and shear stresses 
obtained using the present axisymmetric BEM analysis is compared with 
Reisner's plate theory and a boundary element solution based on 
Reisner's plate theory (Van der Weeen, 1982). General agreement is 
exhibited by the various analyses. The BEM solution of the present 
analysis, deviates slightly from the plate theory results under the 
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origin and near the support. However, this is expected, since plate 
theory is invalid near such conditions. 

The same problem was analyzed with a single region mesh. Results 
were unchanged from those shown. 

6.2.2 Conical Water Tank 

A BEM analysis of a water tank shown in Figure 6.4 is carried out 
using the axisymmetric gravitational particular integrals presented in 
Chapter 3. The hoop stress on the inside surface of the tank 
subjected to hydrostatic pressure is shown in Figure 6.4 along with 
the result from Zienkiewicz (1977). The two results are similar 
although an exact comparison is not possible since geometry data and 
material constants are not given for the later case. In Figure 6.5, a 
comparison of BEM solutions are presented for the hoop stress, on both 
the inner and outer surface (along A-B) of the tank, for two different 
loadings. The first loading is hydrostatic pressure only, and the 
second is a combination of the hydrostatic pressure and self-weight. 
Note that in the lower section the hoop stress at the outer surface is 
compressive, i.e. the bending effect is dominant there. 

One hundred isoparametric quadratic boundary elements are used to 
model the body. The tank is assumed to be constructed of concrete 
with p = 4.65 lbm/ft. 3 , E = 3.472xl0 6 psi, and v = 0.18. The 
thickness of the uniform section is taken to be 1 ft. The time 
required for integration is reduced by dividing the tank into six 
subregions. Two or more elements are used at every interface to 
insure an accurate representation of the variables across the 
interface. The computation time for a single region mesh is one and 
one-half that of the six subregion model, although the computed 
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results are almost identical. A solution to within 5% of the one 
shown is achievable with half the number of elements. 

6.2.3 Rotating Hub 

The axisymmetric BEM analysis of a rotating disk is carried out 
using the centrifugal particular integral discussed in Chapter 3. 
Fourteen boundary elements with twenty-nine nodes are used to model 
the axisymmetric hub (see Figure 6.6). In addition to the boundary 
nodes, stresses at eight internal points are determined. In Figures 
6.6 and 6.7, the contours of tangential stress and equivalent stress 
are shown. Inadequate data is given in Zienkiewicz (1977) for a 
detailed comparison of the similar problem, however, the similarities 
are obvious. 

6.3 INELASTIC ANALYSIS 

6.3.1 Thick Cylinder of Two Materials (with Strain Hardening) 

When a body consists of two or more materials a multi-region 
capability becomes a necessity. The axisymmetric representation of a 
thick cylinder of two materials is shown in Figure 6.8 and the 
variable strain hardening curve for each material is described in 
Table 6.1 below. The load-displacement behavior of the inner and 
outer surfaces is shown in Figure 6.9 for both the direct and 
iterative BEM algorithms based on volume integration. Once again 
excellent agreement is achieved between the two methods with the 
direct method producing a solution up to collapse whereas the 
iterative method does not. The load-displacement curve is initially 
straight while both materials are elastic and starts to bend as the 
inner material begins to yield. At this early stage at which only the 
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inner cylinder has yielded, the evidence of the nonlinearities are not 
appreciable. It is not until the outer material begins to yield that 
substantial nonlinearities develop. 

The cpu times were essentially the same for both methods. 


TABLE 6.1 

Material 1 Material 2 


Yield Stress 

= 10,000 psi 

Yield Stress 

= 36,000 psi 


Plastic 


‘ Plastic * 

Stress (psi) 

Strain 

Stress (psi) 

Strain 

10,000 

0.00 

36,000 

0.00 

30,000 

0.02 

60,000 

0.04 

40,000 

0.04 

60,000 

1.00 

40,000 

1.00 




6.3.2 Steel Pressure Vessel 

An axisymmetrvic, elastoplastic analysis of a vessel subjected to 
internal pressure is shown in figure 6.10. The vessel, constructed of 
steel, has a modulus of elasticity, Poisson's ratio, and a yield 
stress of E = 29.12 x 10 6 psi, v = 0.3, and ' o Q = 40,540 psi, 

respectively. The Von Mises criterion is assumed with no strain 
hardening. Six regions, ninety-nine quadratic boundary elements, and 
twelve quadratic volume cells are used to model the body (Figure 
6.11). Using engineering intuition, cells must be placed in areas 
where yielding is anticipated. The weld connection between the 
spherical shell and branch is of prime concern and it is in this 
region that the cells are situated. The other five regions are 
assumed to remain elastic, and at the end of the analysis it must be 
verified that the stresses in these regions have not violated the 
yield condition. 
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A plot of the vertical deflection of point A with increasing 
pressure obtained from both the variable stiffness and iterative 
plasticity algorithm by volume integration is shown in Figure 6.12 
along with results obtained by Zienkiewicz (1977) using the finite 
element method and test results by Dinno and Gill. The BEM results 
are in excellent agreement with the results obtained by FEM. The 
numerical solutions slightly deviate from the experimental results. 
This variation is due to the idealizations in the numerical analysis, 
such as ideal plasticity and the adoption of Von Mises criterion. 
More importantly, the body is assumed homogeneous when in reality the 
weld and the surrounding region will exhibit greater stiffness. 

It should be noted that the variable stiffness method proceeds 
farther along the curve than the iterative procedure which does not 
converge at these load levels. The cpu times of the two BEM analysis 
were about equal, although on a virtual memory computer (HP9000) the 
real time of the direct plasticity method was greater than the 
iterative procedure due to excessive page faulting. The use of a 
multiregion system in both BEM analysis reduces the computational time 
dramatically, since each nodal equation need only be integrated over 
the surface in which it is contained. Furthermore, the volume 
integration need only be performed over the region containing cells. 

In the past it has been said that BEM analysis should not be used 
on bodies of narrow cross-section, but through sophisticated numerical 
integration and the utilization of a fine mesh, excellent results are 
obtained. 
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6.3.3 Residual Stresses in a Cylindrical Rod 

An axisymmetrie, thermoplastic problem of practical interest 
shown in figure 6.13, regards the residual stresses in a long 
cylindrical rod induced by cooling. In this problem the temperature 
of a rod constructed of 1060 steel is raised gradually to 1250°F and 
then quenched in a brine spray, quickly lowering the temperature to 
80°F. Experimental results obtained by Carman and Hess at Frankford 
Arsenal are given in Boley and Weiner (1960) for the residual stresses 
through the cross-section on the midplane of the rod. The Von Mises 
yield criterion is assumed with a yield value that varies with 
temperature 

-14.3 t + 48710.0 for t < 400 

-18.7 t + 50470.0 for 400 < t <. 775 

-46.3 t + 71900.0 for t > 775 

where or Q = yield stress (psi) and t is current temperature (F°). 
An assumed Biot number of 10 is used in determining the rate in which 
heat is diffused through the rod. 

Utilizing the iterative plasticity algorithm, ten quadratic 
boundary elements and twelve quadratic volume cells were enough to 
produce excellent results for this rather complex problem. Symmetry 
was utilized by the introduction of a roller boundary condition on the 
horizontal bottom face. The residual stresses along this face are 
given in figures 6.14, 6.15 and 6.16. 

6.3.4 Flexible Circular Footing 

This problem concerns a flexible circular footing on an elastic- 
ideally plastic half space which has a modulus of elasticity of E = 
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10,000 kN/m 2 , Poisson's ratio of v = 0.25 and a yield stress of c Q = 
173.2 kN/m 2 (Cu = a Q / -J 3 ) . A Von Mises criterion and associated fiow 
rule is assumed. 

For this axisymmetric analysis, two modeling regions are used. 
One region encloses the anticipated plastic zone and the other defines 
the remainder of the half space which is assumed to remain elastic. 
In Figure 6.17, elements modeling the infinite half space are shown 
starting at the axis of symmetry and continuing along the surface up 
to a finite distance where it is assumed that additional elements 
(modeling the infinite boundary) do not affect the functions in the 
area of interest. (In an argument similar to the three-dimensional 
proof described by Watson (1979), this assumption can be shown to hold 
valid for axisymmetry.) When calculating the coefficients for the 
singular node of the F (traction) kernel via the 'inflation mode 
technique', the modeled region must be completely bounded by a 
surface. Keeping in mind that these coefficients are dependent only 
on the geometry of their respective boundary element, the open region 
can be artificially closed with an arbitrary surface. This arbitrary 
surface is defined with so called 'enclosing elements' and their only 
role is in the calculation of the singular boundary coefficients. It 
is to be understood that the nodal points of the enclosing elements do 
not become part of the boundary system. 

In order to find the convergence of the correct solution, a 
number of cell patterns were constructed for the plastic region. Two 
of these meshes that were fine enough to produce satisfactory results 
are shown in Figure 6.18a and 6.18b. The load-displacement behavior 
at the center of the footing is shown in Figure 6.19 along side the 
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results obtained by Cathie and Banerjee (1980) where the exact 
analytical collapse load for this problem is 6C u . The results of 
Cathie and Banerjee (1980) are stiffer than the results of the present 
paper. The reason for this is the original analysis utilized cells 
with constant variations of initial stress rates in contrast to the 
quadratic shape functions used in the present analysis. More 
importantly, the first investigators calculated stress rates via a 
linear shape function representation for displacement rates (as in the 
finite element method), whereas the present analysis uses an accurate 
integral representation for this. In Figures 6.20 and 6.21 the 
vertical stress and horizontal stress along the axis of symmetry is 
shown. 

It is important to point out a complication associated with the 
idealization of the flexible footing. Any applied loading 
discontinuity can be accurately represented by the boundary element 
discretization since tractions are described on a per element basis. 
However, the boundary stress rates do not as easily admit to a 
discontinuity. In the present analysis the stress rates are 
calculated individually at the node on each adjoining elements using 
the boundary stress calculation and the average value is assumed. 
Therefore, the stress rates are smoothed out across the discontinuity. 
Alternative ways to handle this difficulty are to model the edge of 
the footing as a ramp function across a single element and thus 
avoiding the discontinuity, or a double node can be introduced at this 
point (one for each adjoining cell) and stress averaging eliminated. 
In any case the element(s) containing the troubled point should be 
kept small to better approximate the stress across the discontinuity. 
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6.3.5 Three-dimensional Analysis of a Notch Plate 
The plastic deformation of a notch plate subjected to tension is 
analyzed under plane stress conditions. Four combinations of analyses 
are considered: 

Particular Integral - Iterative 
Particular Integral - Variable Stiffness 
Volume Integral - Iterative 
Volume Integral - Variable Stiffness 

The material properties for the plate are: 

E = 7000 kg/mm 2 
v = 0.2 

<? 0 = 24.3 kg/mm^ (Von Mises yield criterion) 
h = 0.0 

A 90° notch is cut out of the sides of the plate. The maximum to 
minimum width ratio is 2 and the thickness is 6/10 of the maximum 
width. A quarter of the plate is discretized in two subregions as 
shown in figure 6.22. The region, containing the notch, has thirty 
quadratic boundary elements. The inset in this figure shows six 
(twenty-node) isoparametric cells (with sixty-eight distinct cell 
nodes) which are used in the volume integral based analysis. In the 
particular integral analysis, particular integrals are defined using 
the boundary nodes and three interior nodes (corresponding to the mid- 
side nodes of the cells) for a total of ninety-two particular integral 
nodes. The second region has sixteen quadratic boundary elements. 
Boundary conditions on the front and back faces are assumed traction 



In figure 6.23, the stress-strain response on the mid-plane of 
the root is given for the four (3-D) BEM analyses and compared with 
the two-dimensional plane stress and plane strain BEM solutions 
obtained by Raveendra (1984). The three-dimensional results are in 
good agreement with one another, and fall between the two-dimensional 
solutions, closer to the plane stress result, as one would expect. 
For solutions above the load level 2 <J m /a 0 = 1.0 requires a finer mesh 

around the notch, particularly in the volume cell discretization. 

‘ * 

Figure 6.24, shows another notch plate mesh for a particular 
integral analysis. The boundary discretization is the same, but ten 
additional particular integral nodes are added in the interior (101 
particular integral nodes, total). The bottom face becomes a plane 
of symmetry by applying a roller boundary condition. In order to keep 
the dimension of *the plate proportional to the first analysis, the 
thickness of the mesh must be reduced by one-half, since symmetry is 
assumed by virtue of the additional roller boundary condition. The 
stress-strain response at the root, obtained for this mesh using the 
variable stiffness algorithm, is shown in figure 6.25. Results for 
the three nodes across the thickness of the root are compared with the 
two-dimensional solution (Raveendra, 1984). Once again, the three- 
dimensional results lie between the plane stress and plane strain 
solution. 

6.3.6 Three-dimensional Analysis of a Perforated Plate 

The plastic deformation of a perforated plate in tension is 
analyzed under plane stress conditions. The volume distribution of 
initial stress is represented using either twenty-noded isoparametric 
volume cells or the point based particular integral. Each 
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representation is coupled with either the iterative or the variable 
stiffness solution process, leading to a total of four distinct 
algorithms. The material properties for the plate are: 

E = 7000 kg/mm 2 
v = 0.2 

or Q = 24.3 kg/mm 2 (Von Mises yield criterion) 
h = 224.0 kg/mm 2 

The diameter of the circular hole, at the center of the plate, is one- 
half the width, and the thickness is one-fifth the width. A quarter 
of the plate is discretized in two subregions, as shown in figure 
6.26. The first region, containing the root of the plate, has thirty 
quadratic boundary elements. For the volume integral based analysis 
the domain of this region is discretized using nine (twenty-node) 
isoparametric cells. In the particular integral analysis, particular 
integrals are defined at points corresponding to the cell nodes used 
in the volume integral analysis. The second region has twenty-three 
boundary elements. No volume discretization or definition of 
particular integrals is required in this region since it remains 
elastic throughout the analysis. Boundary conditions on both the 
front and back faces are assumed traction free. 

This problem was previously analyzed experimentally by Theocaris 
and Marketos (1964) and by Zienkiewicz (1977) using the finite element 
method. The results obtained by the boundary element analysis is 
compared to these results in figures 6.27 and 6.28. The stress-strain 
response at the root of the plate is shown in figure 6.27. The 
results obtained using the various BEM algorithms show good agreement 
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with one another and with the variable stiffness FEM analysis. 
Differences between the iterative and variable stiffness BEM 
formulations are much less significant than the difference in the two 
FEM algorithms. The volume integral base BEM procedure exhibits 
greater stiffness than the particular integral method. In figure 
6.28, the stress distribution between the root and the free surface of 
the specimen is shown for a load of 2 ® m /o 0 = 0.91. Once again, 
excellent agreement is obtained among the four boundary element 
analysis. In order to evaluate the degree of convergence of the 
results, a mesh with sixteen volume cells (or the particular integral 
equivalent) was studied. The results were unchanged from those shown. 

The present analysis was carried out on the Cray-1 computer. The 
CPU times for the four algorithms were: 


Particular Integral/Iterative 
Particular Integral/Variable Stiffness 
Volume Integral/Iterative 
Volume Integral /Variable Stiffness 


= 254 sec. 
= 358 sec. 
= 272 sec. 
= 361 sec. 


6.3.7 Two-dimensional Analysis of a Perforated Plate 

The perforated plate of the previous problem is analyzed using 
the particular integral based, inelastic, two-dimensional, plane 
stress analysis. Three discretizations, shown in figure 6.29, are 
used to model the plate. Each mesh is divided into two subregions. 
The initial stress distribution in the inelastic region is defined 
using the particular integral representation. Twenty, forty-one and 
sixty-five nodes are used, respectively, in these discretizations to 
define the particular integrals. 
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A solution is obtained using the variable stiffness method, and 
the results are compared with the particular integral based, variable 
stiffness results of the previous three-dimensional analysis. The 
stress-strain response at the root of the plate is shown in figure 
6.30. The results obtained for the two-dimensional discretizations 
are in good agreement with the three-dimensional solution. The 
results of the more refined (65 particular integral nodes) mesh, 
however, exhibits a smoother response and follows the three- 
dimensional solution more closely than the other two-dimensional 
results. The axial stress distribution across the plate, from the 
root to the free edge, is shown in figure 6.31. All two-dimensional 
results vary from the three-dimension solution, however, overall 
agreement is observed. 

6.4 CONCLUDING REMARKS 

In this chapter, a range of problems of practical interest have 
been analyzed using the axisymmetric, two- and three-dimensional, 
thermal, elastic and inelastic boundary element formulations of 
previous chapters. In addition to displacement and traction boundary 
loads, body forces due to self-weight, centrifugal, and thermal loads 
are considered. The solutions obtained in the analyses were compared 
to known solutions when they exist, and good agreement was found in 
most cases. Many of the problems were analyzed using various 
alternative procedures to illustrate the relative accuracy of these 
methods. Also demonstrated was the multi-region capability of the 
computer code which enables one to model a body in substructured 
parts; thus dramatically reducing the cost of the analysis. 
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Figure 6.1a 

Clamped Circular Plate 



Figure 6.1b 

Axisymmetric Mesh of a Circular Plate 
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Figure 6.2 

Radial and Tangential Stress through a Circular Plate 
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Figure 6.3 

Shear Stress though a Circular Plate 
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Figure 6.4 

Conical Water Tank subjected to Hydrostatic Pressure 
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Figure 6.5 

Hoop Stress in a Conical Tank due to Hydrostatic Pressure and Self-weight 
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Figure 0.6 

Contours of Tangential Stress in a Rotating Hub 
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Figure 6.8 

Axisymmetric Mesh of a Thick Cylinder of Two Materials 
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Figure 6.9 

Radial Displacement of a Two-Material Thick Cylinder under 
Internal Pressure (Plane Strain) 





Figure 6.11 

Discretization of Pressure Vessel Mesh by Region 
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Figure 0.12 

Vertical Deflection of Point A of the Pressure Vessel with Increasing Pressure 
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Figure 6. IS 

Axisymmetric Mesh of a Cylindrical Rod 
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Figure 0.14 

Residual Radial Stress in a Cylindrical Rod 
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Figure 6.15 

I Residual Tangential Stress in a Cylindrical Rod 
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Figure 6.16 

Residual Axial Stress in a Cylindrical Rod 
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Figure 6.17 

Boundary Discretization of an Axisymmetric Halfspace under a Flexible Circular Footing 
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Figure 6 . 18 a 

Mesh 1: Volume Discretization of the Plastic 
Region under a Circular Footing 



Figure 6 . 18 b 

Mesh 2: Refine Mesh of the Plastic Region 
under a Circular Footing 
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Figure 6.19 

Load-Displacement Behavior under a Circular 
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Figure 0.20 

Vertical Stress at the Origin under a Circular Footing 


DEPTH / RADIUS OF FOOTING 



.75 .00 -.75 -1.50 -2.25 -3.00 -3.75 -4.50' -5.25 -6:00 -S.75 

RADIAL STRESS / Cu* 


Figure 6.21 

Radial Stress at the Origin under a Circular Footing 
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Figure 6.22 

Boundary and Volume Discretization of a Three-dimensional Notch Plate 
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Figure 6.23 

Stress-Strain Response (on the Mid-plane) at the Root of a 
Three-dimensional Notch Plate 
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Figure 6.24 

Three-dimensional Discretization of a Notch Plate for Particular Integral Analysis 



Figure 6.25 

Stress-Strain Response at the Root of a Notch Plate, 
Particular Integral - Variable Stiffness (Plane Stress) Analysis 
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Figure 0.20 

Three-Dimensional Mesh of a Perforated Plate 
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Figure 0.27 

Stress-Strain Response at the Root of a 3-D Perforated Plate 
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Figure 6.28 

Stress Distribution across a 3-D Perforated Plate near Collapse Load ( a ^ 2cr m /a 0 = 0.91) 
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Figure 0*30 

Stress-Strain Response at the Root of a Perforated Plate (Plane Stress) 
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Figure 6.31 

Stress-Distribution across a Perforated Plate (Plane Stress) 
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CHAPTER SEVEN 


CONCLUSIONS AND RECOMMENDATIONS 
7.1 GENERAL CONCLUSIONS 

Considerable effort was aimed at extending all formulations 
presented in this dissertation to their axisymmetric form. * The 
inelastic axisymmetric BEM analysis of this work represents the most 
advanced implementation of its kind. The initial stress expansion 
technique was particularly helpful in axisymmetry to overcome the 
difficulty in calculating the coefficient of a singular point near the 
origin. 

A direct, variable stiffness type inelastic solution algorithm 
was, for the first time, implemented in a general purpose, multiregion 
computer code. Although the computational time of analysis for this 
method is greater than that of the conventional iterative procedure, a 
distinct advantage of the variable stiffness method is its ability to 
produce a solution close to the collapse state of stress. 

New boundary element formulations, based on particular integrals, 
were introduced for the treatment of body forces and nonlinear 
effects. The method eliminated the need for volume integrals or extra 
surface integrals to account for these effects. Furthermore, the 
method is applicable to other BEM analysis which involve an 
inhomogeneous differential equation. 

Finally, inhomogeneous formulations were presented for elastic 
and (for the first time) inelastic media. The formulations accounted 
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for thermally induced inhomogeneities and those resulting from natural 
variation in material parameters. The equation system of this 
analysis is arranged in the same form as the homogeneous-material 
system. This allowed existing inelastic solution algorithms to be 
employed without modification. 

7.2 RECOMMENDATIONS FOR FUTURE WORE 

In order to facilitate future research based on the findings of 
the present work, the following recommendations are presented: 

1 . The axisymmetric analysis of the present work is capable of 
handling axisymmetric (boundary and body force) loading 
only. Although the axisymmetric elastic analysis for 
arbitrary loading exists, to date, the inelastic formulation 
does not, and therefore, should be developed. 

2. Further refinement for efficiency is needed in the variable 
stiffness inelastic solution algorithm. At present, a 
solution of a new system matrix is required at every load 
step. A new strategy, utilizing a re-solution of the 
previous increment's decomposed matrix system at alternate 
load steps, may be possible. 

3. The variable stiffness method is more stable than the 
iterative algorithm near the collapse state of stress, 
however, the iterative method is capable of plastic strain 
softening where the present variable stiffness 
implementation is not. It is conceivable that a hybrid 
method could be developed that would incorporate the best 
features of both methods. 
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4. The present work has been more concerned with the 
development of BEM as an instrument for solving inelastic 
problems than with the development of inelastic models for 
BEM. This was obvious since only the Von Mises yield 
criterion (with hardening) was employed. Now with the 
present inelastic BEM analysis refined to such an advanced 
stage, more realistic material models need to be included 
for the study of more complex materials. 

5. The analyses in this dissertation are limited to problems of 
material nonlinearities. Problems such as metal forming 
require the modeling of geometric non linearities associated 
with finite displacement theory. BEM formulations should be 
developed based on this theory. An extended form of the 
particular integral based method can be derived for this 
purpose or the conventional volume integral approach can be 
used. 

6. Although the present functional form of the global shape 
function, chosen for the approximation of the particular 
integrals, produced satisfactory results, further research 
still should be devoted to the discovery of new functions. 

7. The axisymmetric and two-dimensional inhomogeneous 
formulations developed in this dissertation should be 
extended to three-dimensional analysis. However, the 
extension to three-dimensions greatly increases the overhead 
of the analysis, and therefore, an efficiency study should 
be carried out to determine its feasibility. 
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8. The variable stiffness algorithm, the particular integral 
formulation, and the inhomogeneity formulation all have one 
thing in common. All involve a large amount of matrix 
manipulations. Therefore, to maximize the efficiency of 
these new methods, the computer code should be vectorized to 
exploit the state-of-the-art computer technology. 
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APPENDICES 


APPENDIX I - NUMERICAL INTEGRATION 
I. A Shape Functions 

When a body is discretized, as described in Chapter 2, shape 
functions are used to describe the geometry and field variables across 
the boundary element (or cell) in terms of nodal variables. Geometry 
coordinates; displacements, tractions, and initial stresses are 
expressed in terms of shape functions as 

x i = N®(£) x® 

u ± = N®( £) u® 

t i = N®( £) t® 

where 

N® and represent the shape function, 
a and 0 the order of the shape function, 

5, and 1^2 are intrinsic (local) coordinates, and 
the bar indicates the nodal values. 

Some common shape functions used in this dissertation are given 
below. Additional shape functions can be found in text books, (e.g. 
Banerjee and Butterfield, 1981). 


220 



Three-noded, one-dimensional shape functions (Figure A.l): 

M 1 (n) = | nd-Hi) 

m 2 (t,) = - \ nd-n) 

M 3 (t,) = (l+Tl)(l- n ) 

Eight-noded , two-dimensional shape function (Figure A. 2): 

M = ^(1+s^Tij) (l+s^Hj) ( s 0 itli+s a 2ti2 - l) = 1*2, 3, 4 

M a (n 1 ) = |(l+S a 2Tl 2 Hl-^i) a = 5,7 

M a (ti 1 ) = |d + s ol Tl 1 )(l-T l ^) a = 6,8 

where s ^ takes the sign of coordinates. The shape functions for 
six-noded triangular elements can be obtained by collapsing the 
quadrilateral to triangle. 

Twenty-noded , three-dimensional shape function (Figure A.3): 

= I (1+S ol^l )(1+s a2 T '2 )(1+s a3 T '3 )(s al 1 ^l +s a2 T l2 +s a3 T »3 _2) 

a=l . . 8 

M ° (n l ) = 4 -( 1 " t, l )( 1 +s o2^2 )( 1 +s a3 1 l3 ) a=9..12 

M a (n 1 ) - |(l-n 2 ) < 1 +s ol T 'l ) < 1 + 8 a 3 T > 3 ) 0 -IS..I 6 

M ° ( n 1 ) = |( 1 -^3)(l + s ol n 1 )(l+s a 2 Ti2) c= 17 . .20 

where s q1 assumes the sign of coordinates. 
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I.B Jacobian Transformations 
Boundary element mapping: 

A curvilinear (2-D) line element can be mapped into one- 
dimensional space using the following Jacobian transformation 


dC(x) = J L (n)dti 


where 


a Xi a Xi 1/2 

= L FT ail J 


i = 1,2 for a two-dimensional 
curvilinear element 


Similarly, a curvilinear (3-D) surface element can be mapped into two- 
dimensional space using the following Jacobian transformations 


dS(x) = J s (i) 1 ,Ti 2 )dTi 1 dT |2 


where 




r .2 . .2 . 2 , 1/2 

[dj. + d 2 + d 3 ] 


d, = < 


9X 2 

8x 3 _ 

3x 2 

3X 3 


dx] 2 

3ti 2 

9^2 


3x 3 

3 Xl 

3x 3 

Sil^ 

3t1 2 

dt l 2 


3Xi 

3X2 


3x 2 

3^2 

3n 2 

3n 2 

9h 2 


Volume cell mapping: 

Two- and three-dimensional volume cells can be mapped into two- 
and three-dimensional space, respectively, using the following 
Jacobian transformations 
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dA(x) - ^1^2 

dV(x) = J(n 1 ,n 2 *ti 3 ) dn^dri 2 T l 3 


where 


J(\) 


3x. 

Ila“ll 

8,, j 


i,j = 1,2 for two-dimensions 
i,j = 1,2,3 for three-dimensions 


I.C Gauss-Legendre Formula 

‘ • 

In one dimension, the application of. Gauss-Legendre formula 
(Stroud and Secrest, 1966) is expressed as 


A A 

) f(x)dx = 5 w a f(x a ) 

-1 a=l 

where 

f(x) is the prescribed function, 

w a is the weight at sampling point a, 

s a is the abscissae of the sampling point, and 

A is the order of the integration rule. 


Two- and three-dimensional integrals are evaluated through 
repeated application of the above formula: 

f l 1 A B 

J f f (x,y)dxdy = \ \ f(x a ,x b ) 

-1-1 a=l b=l 


and 


J { If (x,y,z)dxdydz \ w a w b w c f (x a ,y b ,z c ) 

^ —A — 1 wr-« 


ABC 

\x l 

a=l b=l c=l 


223 



APPENDIX II - TWO- AND THREE-DIMENSIONAL KERNEL FUNCTIONS 

The kernel functions for the two- and three-dimensional boundary 
integral equations for displacement and stress are presented below. 
In these definitions the parameter d=2 is used for two-dimensions 
(plane strain) and d=3 for three-dimensions. The indicial notation 
ranges from 1 to 2 for two— dimensions , and from 1 to 3 for three- 
dimensions. The two-dimensional plane stress kernel functions are 
derived from the plain strain function by using a modified material 
parameter v = /(1+ v). 


II. A Displacement Equation 

Uj(lj) = f s [G i j(x,?)t i (x) - F^(x,5)u i (x)]ds(x) 

+ f y [G i j(x,5)f i (x) + B ik j(x,5)c° k (x)]dv(x) 

in which G . is the fundamental point force solution due to Kelvin 

i J 

(Love, 1944) 

0 

G ij (x,?) = — [ c 2 6 ij £ (d-2) + (d-3)ln(r)] + 


where 


_ 1 

C 1 8np(l-v) 
C 2 = (3-4 v) 


y i " x i"*i 

r 2 = y^i 
Z i - y i /r 

and d is the dimensionality of the problem. 
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G^j represents the displacement u^(x) at * in direction i due to a 
point force e^) a t 5 in direction j, 

i.e. u ± ( x ) = 

It should be noted that the two-dimensional solution includes an 
arbitrary constant term in addition to the expressions provided by 
the above equations since# unlike the three-dimensional case where G^ 
vanishes as r-»®, the expression for G^ does not vanish as r-*® for 
two-dimensional problems. 

The kernel represents the strains corresponding to the 

fundamental displacement solution and is given as 

B ij k <x.5> = - (d-l) [ c 4 ( 8 jk z i + 8 ik z j ) “ 8 ij z k + dz i z j 2 kl 

id-ur 


where 


c 4 = ( 1-2 v) 

The kernel represents the tractions on an incline normal 
n^(x) due to a point force: 


F ij(x,5) = 


(d-l)r 


WT) [ C 4 (n j z i" n i z j ) + (C 4 8 ij + dz i z j )z m n m] 


II. B Stress Equation 

<*^(5) = / s [G^jkCx.Otj^Cx) - F^j k (x,?)Ui(x)]ds(x) 

+ f v tG ijk (x '5>fi(x) + Bj pjk (x.$)®ip]dv 

where 

c 3 

G iJ k (x»5> = - (d-l) [ c 4 ( 8 ij z k + 8 ik 2 j " 8 jk z i ) + dz i z j z k] 

id-ijr 
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and 


F ijk (x '^ ” ( d _ 1)r d ln i (C 6 6 jk " dC 4 z j z k > " n j (C 4 8 ik + d(v)z i z k ) 

" V C 4 8 ij + d ^> x i 2 j> - dz m n m IC 4 8 jk z i + ( v )8 ik z j 
' (v) 8 ij z k " U + »*i*j* k n 

C 3 

B ipjk (z ' 5> = (dl)r d tC 4 (8 ip 8 Jk “ 8 ij 8 kp " 8 ik 8 jp “ d8 jk z i z p ) 

" d8 ip z j z k - d(v)(8 ij z k z p + 8 jp z i z k + 8 kp z i z j + 8 ik z j z p> 

+ d(d+2)z i ZjZ k z p ) 

= -1 

C 3 4n ( 1- v) 
c = — -=M 

L 5 2n(l- v) 

C g = (1-4 v) 

The integration of the kernel is strongly singular and must 

be treated as a Lebesgue integral. This term can be decomposed into 
Cauchy principle-value integral and a free term 

/„ Bj jkl ( x ,5)dv(x) - /„ 5j Jkl (x,?)dv(x) + Jj jkl (?> 

where 

C 

J Ipjk - d(dl2) - <-Md 2 -4)]6i j6kp + [l-(v,(d+l)8 lp Sj k ]] 

c - -=l— 

0 8 (1-v) 
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APPENDIX III - AXISIHHETRIC KERNEL FUNCTIONS 


Here, the displacement and stress kernel functions used in the 
axisymmetric boundary element analysis of bodies subjected to 
axisymmetric boundary loads and initial stress body forces are given. 

The kernel of equation (III.2) is the ring source solution, 

first derived by Kermanidis (1975) and Cruse et. al. (1977). Other 

authors Baker and Fenner (1983) and Tan (1984) have defined as 

the transpose of the ring solution. In the following r and z 

correspond to the field point, and R and Z are the integration (Gauss) 

3 8 

points. Note, although ~ jz holds true for axial (z) direction 

in axisymmetric kernels, such is not the case for the radial (r) 
direction, i.e. “ £ - |p. 

A 2jiR term, which appears after integration in the O 
direction, has been absorbed in the kernelv therefore, the numerical 
integration is performed over a curve dC . Some other authors Bakr 
and Fenner (1983) and Cruse et. al. (1977) do not absorb this term in 
their kernels. 


III. A Displacement Equations 


u r (r,z) 

u z (r,z) 


■ J< 

C 


r r zr 


G rz G zz 


F F 
r rr r zr 


F rz F zz 




|U„ 


> dC(R,Z) 


* J 

B 11 

B 21 

B 31 

B 41 

A 

B 12 

B 22 

B 32 

B 42 


J rr 

o 

J zz 

■f. 


dA(R,Z) 


(III.l) 
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in which 


G ij = A ij K + B ij E = r ' 2 (III. 2) 

being the displacement vector due to a ring load 

intensity e^, i.e., u^x) = G ij (x 'S )e j 
K = K(m) is the complete Elliptic Integral of the first kind, 
and 

E = E(m) is the complete Elliptic Integral of the second kind. 


F . 
n 


r 8t >i ,°j ± , 

L C 1 3R + c 2 ( R 3Z J n r + 11 l 3Z 8R J n 2 


zi 


8G , G . 3G . 8G . 6G , 

r zi . , ri . ri,i r n . zii 

L C 1 ~YT + c 2 ( ~ + TT ) J n z + » ITT + IFJ n r 


i=r,z 


aG 


rr 


9G 


zr 


3 11 “ dR C 21 “ ez B 31 az T 3R 


aG rr + aG zr 


„ _ rr 

B 41 “ R 


aG 


rz 


3G 


12 aR B 22 az D 32 az ' aR 


zz 


B, 


3G rz + aG zz 


_ zjrz 
B 42 ~ R 


where 


, 3A • j 3B, j atf 

1 a ir + — J p . » — + d 

3 R - 3 R K aR E + A ij 3 R B ij 9 R 


9E 


90. . a A 


3B • 


il = !^l K + ZM E + 

ry *rr * rj £j“T 


az az 


az 


0K 3E 

E + A ij az + B ij az 


i,j = r,z 


n = n (K,Z) is the normal in r direction 
“r r 

n z = n z (R,Z) is the normal in z direction 
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3A Z 

Z2 

3Z r 2 a zz 


SB 4«Z 

— = — = -r- - d a B 

.2,, 3 ZZ 


3Z 


pH 


H = [(R + r) 2 + l 2 ] 1 ^ 2 X,(i = Lame constants 


Z - (Z - z) 


v = Poisson's ratio 


P 2 = [(R - r) 2 + l 2 ] 


r 
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= field points 
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8np(l-v) 
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= integration (Gauss) points 
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d 5 = N + 2l 2 


— 2 

S 1 = c 3 (r + r) + ~~(H 
P 4 


2R 2 ) 


s« = 


c 3 " 


2Z 2 Rr M_ 

n 4 2 
P P 


III.B Stress Equation 
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III.C AxisyametPic Jump Term Tensor, ^*^ 1 ( 5 ) 

The jump term (or free term) that is associated with the Lebesgue 
domain integral of the stress equation is given below. This jump term 


corresponds to the axisymmetric 'initial stress' formulation defined 
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III.D Elliptic Integrals 

The complete elliptic integral of the first (K(m)) and second 

(E(m)) kind with modulus m can be defined as 
n/ 2 

K(m) = 


I 


dO 


[1-rn sin 2 0] 1/2 


n/2 


E(m) = J [l-m sin 2 ©] 1/2 d0 

The modulus corresponding to the axisymmetric kernel is given as 


m = — t— m. = l - m 
n 1 

For numerical calculations, the integrals are approximated using the 
following polynominal approximation (Abramowit 2 and Stegun, 1974) for 
(0 < m < 1) . 
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K = K(m) = £ [ a i m i + + 
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E = E(m) = 1 + ^ [ c i m i + d i m i 
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|e(m)l is <. 2 x 10 ® and a^, b^, c i and d A are constants defined 

below: 
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Figure A.l 

Two-dimensional Boundary Element 



Figure A.2 

Three-dimensional Boundary Element 
(or Two-dimensional Volume Cell) 



Figure A.3 

Three-dimensional Volume Cell 
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